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1  Summary 

This  report  summarizes  findings  from  a  three-year  research  effort  aimed  at  enhancing  the  understanding  of 
primary  atomization  processes  and  their  role  in  liquid  rocket  engine  combustion  instabilities.  A  new  fully 
nonlinear  primary  atomization  model  has  been  developed  to  provide  fundamental  information  for  droplet 
sizes  for  arbitrary  injection  conditions.  This  model  makes  use  of  no  empirical  parameters  thereby  providing 
a  general  tool  capable  of  addressing  the  influence  of  various  design  changes  on  the  primary  atomization 
process.  The  model  shows  good  agreement  with  limited  experimental  data  as  the  experiments  required  to 
confirm  the  results  are  very  challenging. 

We  have  also  developed  a  homogeneous  model  for  addressing  two-phase  flow  inside  submerged  regions  of 
coaxial  injectors.  This  unsteady,  3-D,  Navier-Stokes  calculation  has  been  used  to  simulate  flows  in  the  Space 
Shuttle  Main  Engine  injector  successfully;  although  limited  data  are  available  results  show  a  good  correlation 
with  fundamental  flame  flickering  frequencies  in  recent  experiments  at  DLR.  The  model  has  been  used  to 
characterize  the  amplitude  and  frequency  of  the  instability  over  a  wide  range  of  operating  conditions. 


2  Research  Objectives 

Liquid  rocket  engine  combustion  chambers  represent  the  most  intense  energy  release  (on  a  per  unit  volume 
basis)  of  any  man-made  device.  Unfortunately,  this  situation  makes  the  combustors  particularly  sensitive  to 
instabilities  -  a  problem  that  has  plagued  the  industry  virtually  since  its  inception.  The  objectives  of  this 
research  program  are  to  assess  the  role  of  the  atomization  process  in  contributing  to  instabilities,  to  assess 
the  potential  frequencies  of  the  instabilities  and  thereby  to  provide  mechanisms  to  control  the  onset  of  these 
damaging  phenomenon.  The  following  section  details  the  status  of  the  research  in  achieving  these  goals. 

3  Status  of  the  Research 

3.1  A  Fully  Nonlinear  Model  for  Atomization  of  High-Speed  Jets 

A  nonlinear  model  has  been  developed  to  assess  the  time  evolution  of  an  axisymmetric  liquid  jet  using  a 
boundary-element  method.  Vorticity  transported  from  the  boundary  layer  in  the  orifice  passage  to  the  free 
surface  is  modeled  using  a  potential  ring  vortex  placed  at  the  orifice  exit  plane.  The  vortex  strength  is 
uniquely  determined  from  the  Kutta  condition  and  information  regarding  the  boundary  layer  thickness  at 
the  orifice  exit  plane.  It  is  shown  that  primary  breakup  can  occur  even  in  the  absence  of  the  gas  phase. 
Using  a  secondary  stability  analysis  after  Ponstein  [1],  the  size  of  the  droplets  is  estimated  based  on  the  size 
of  the  ring- type  structures  shed  from  the  periphery  of  the  jet.  Computed  droplet  sizes  are  in  reasonable 
agreement  with  experimental  data,  although  turbulence  effects  obscure  some  comparisons. 

The  complete  simulation  of  the  Hoyt  and  Taylor’s  jet  [2]  is  shown  in  Fig.  1.  The  jet  structure  is  initially 
assumed  to  be  a  simple  cylinder  with  a  hemispherical  tip.  and  its  evolution  is  simulated  via  time  integration. 
The  simulation  is  completed  at  t  —  5.0.  A  slight  ‘swelling’  is  observed  at  t  ~  1.0  and  the  fluctuation  of  the 
jet  surface  is  seen  at  t  >  2.0.  The  velocities  induced  by  the  bound  vortex  are  large  enough  to  penetrate  the 
jet  surface  and  it  results  in  the  primary  atomization.  It  should  be  noted  that  most  liquid  ligaments,  pinching 
from  the  jet  surface,  are  in  the  ‘rollup’  motion  in  counterclockwise  direction  while  the  mean  velocity  of  the 
ligament  is  in  the  streamwise  direction.  The  counterclockwise  rollup  motion  is  a  strong  evidence  that  the 
boundary  layer  instability  is  the  fundamental  cause  of  the  primary  atomization.  The  mean  velocity  of  most 
droplets  are  in  streamwise  direction  as  the  droplets  motion  propagates  along  with  the  main  jet  stream,  the 
most  dominant  convective  source. 

It  is  well  known  that  the  droplet  size  varies  significantly  within  the  atomization  regime.  Wu  et  al  [3] 
reported  the  droplet  size  variation  with  U  for  the  turbulent  water  jet  into  air.  Hoyt  and  Taylor’s  experiment 
had  been  carried  out  for  the  Bernoulli  pressure  A P  <  60  psi ;  no  result  with  higher  A P  is  reported  [2,  4,  5]. 
However,  we  hypothesized  the  increase  in  A P  for  the  Hoyt  and  Taylor  jet.  The  result  is  taken  the  jet  speed 
up  to  U  —  40  m/s  which  corresponds  to  A P  pa  116  psi  or  slightly  higher  in  order  to  account  for  some 
pressure  loss  within  the  nozzle.  The  final  result  for  the  Hoyt  and  Taylor’s  case  [2]  is  shown  in  Fig.  2.  Using 
the  methodology  employed  in  the  previous  section  with  no  calibration  constants,  the  model  predicts  the 
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Sauter  Mean  Diameter,  SMD,  with  reasonable  accuracy.  As  shown  in  Fig.  2,  there  is  a  steep  gradient  at  a  jet 
speed  around  U  «  20  m/s  for  the  Wu  et  al.  [3]  turbulent  jet  experiment.  Our  model  result  overlaps  with  that 
obtained  by  Hoyt  and  Taylor.  It  is  interesting  to  observe  that  Wu  et  al.’s  data  is  also  similar  to  our  result 
and  that  of  Hoyt  and  Taylor  at  about  U  %  20  m/s.  While  Wu  et  al.  noted  this  as  the  region  of  ‘uncertainty’, 
it  is  possible  that  the  rollup  motion  was  competing  with  the  turbulence  and  thus  the  perceptible  effect  of 
the  rollup  motion  appears  as  shown  in  Fig.  2. 

For  U  20  m/s,  differences  between  the  calculations  and  the  experiments  emerge.  It  is  known  that  linear 
analysis  [6]  overpredicts  the  droplet  size  (by  less  than  20%)  because  it  neglects  the  satellite  droplet  mass  due 
to  the  nonlinear  effect  [7,  8]  which  yields  the  multiple  crests  per  wavelength.  Ponstein’s  equation  [1]  is  a  linear 
analysis  and  thus  it  also  tends  to  overpredict  the  droplet  size.  However,  a  20%  or  smaller  SMD  difference 
does  not  explain  the  difference  we  see  in  Fig.  2  at  a  higher  jet  speed.  This  is  due  to  the  fundamental  difference 
between  the  boundary  layer  instability  jet  and  the  turbulent  jet:  the  boundary  layer  instability  jet  is  scaled 
by  the  momentum  thickness  [9]  and  the  turbulent  jet  is  scaled  by  the  Kolmogrov  length  scale,  lk  or/and 
turbulence  eddy  characteristics  length  of  kinetic  energy,  [3].  Wu  et  al.  derived  the  empirical  formula  using 
the  ‘surface  kinetic  energy’  argument  which  gives  SMD  scaled  by  ~  1  /UlA8.  Thus  the  governing  length 
scale  (i.e.  lk  and  /,)  decrease  significantly  at  about  U  ~  20  m/s.  On  the  other  hand,  the  SMD  of  the 
boundary  layer  instability  jet  is  scaled  by  &2 :  SMD  -  l/U0  5  and  thus  its  change  with  respect  to  U  is 
relatively  moderate  as  shown  in  Fig.  2.  In  Fig.  3,  the  jet  becomes  more  stable  (or  less  atomization)  as  U 
increases.  This  is  consistent  with  Rupe’s  observation  [10].  As  U  increases,  the  nozzle  exit  velocity  profile 
becomes  flatter  which  contains  less  circulation.  One  may  imagine  the  limiting  case  of  such  condition:  U  — ►  oo 
(Red  ->  oo)  and  thus  inviscid  flow  with  the  perfectly  flat  velocity  profile  which  contains  no  circulation  at  the 
nozzle  exit.  In  this  limiting  case,  the  jet  is  unconditionally  stable  unless  it  is  perturbed  by  other  instability 
mechanisms.  If  the  circulation  at  the  nozzle  exit  increases,  the  jet  is  more  unstable. 


3.2  Hydrodynamic  Instability  of  Shear  Coaxial  Jet  in  a  Recessed  Region 

In  many  coaxial  injectors,  liquid  flowing  down  a  central  post  is  atomized  by  a  high-velocity  gas  passing 
around  the  outer  annulus.  In  many  applications,  the  liquid  post  is  submerged  somewhat  from  the  orifice 
exit  plane  to  provide  flame  holding  in  combustion  systems  such  as  liquid  rocket  engines.  The  submergence 
provides  for  oscillations  of  the  liquid  jet  within  the  confines  of  the  submerged  region.  These  disturbances 
could  couple  to  the  dynamics  of  the  jet  breakup  process  and  potentially  provide  amplification  of  oscillations 
within  the  combustion  chamber.  Combustion  instabilities  of  this  nature  can  have  severe  impact  on  the 
performance  of  the  engine  and  can  in  some  cases  lead  to  catastrophic  failures. 

In  general,  the  combustion  instability  with  high  frequency  is  categorized  into  acoustic  instability  and 
hybrid  instability.  The  acoustic  instability  shows  dominant  wave-type  oscillation  in  the  main  chamber, 
but  is  independent  of  the  feed  system.  With  the  hybrid  form  of  instability,  the  wave  character  of  the 
oscillation  is  strongly  coupled  between  the  feed  system  and  the  combustion  chamber.  Hutt  and  Rocker  [11] 
also  investigated  the  high  frequency  combustion  instability  associated  with  coaxial  injectors.  They  classified 
the  instability  phenomena  in  the  chamber  as  injection-coupled  and  intrinsic  mechanism.  The  injection 
coupling  implies  chamber  pressure/temperature  variation  as  a  key  contributor  in  the  change  of  flow  dynamics 
through  the  injector.  In  the  other  hand,  the  intrinsic  mechanism  occurs  in  the  flowfield  due  to  its  own  flow 
dynamics  with  negligible  feed  system  effect.  It  should  be  noted  that  injection  coupling  is  never  independent 
of  the  intrinsic  subprocesses,  such  as  atomization,  propellant  heatup,  vaporization,  and  mixing  because  these 
processes  determine  the  relationship  between  the  injector  response  and  the  chamber  response. 

However,  in  hydrodynamics  point  of  view,  for  the  flow  at  "high  velocity,  most  of  researchers  agree  that 
the  principal  source  introducing  instability  to  the  jet  is  from  aerodynamic  forces  arising  from  the  interaction 
of  the  liquid  jet  with  the  surrounding  gas  flow.  Reynolds  and  Weber  numbers  are  generally  very  high  in 
these  atomizers  and  aerodynamic  forces  are  several  orders  of  magnitude  larger  than  capillary  forces.  The 
interaction  between  the  liquid  and  gas  phases  mainly  comes  from  different  velocities  of  each  phase.  The 
velocity  discontinuity  in  a  homogeneous  fluid  results  wave  growth  on  the  interface,  which  is  a  common 
Kelvin- Helmholtz  instability. 

Mayer  and  his  research  group  [12,  13,  14]  have  done  a  significant  amount  of  work  on  the  coaxial  injector 
in  terms  of  the  combustion  instability.  They  claimed  the  initiation  of  the  jet  surface  deformation  was  due  to 
internal  liquid  turbulence  transforming  energy  between  phases  in  forms  of  eddy  structures,  approximately 
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Figure  1:  Jet  evolution  for  conditions  consistent  with  Hoyt  and  Taylors  experiment  [2].  Annular  ligaments 
which  have  pinched  off  from  the  domain  are  not  shown  to  improve  clarity. 

a  size  of  10-30%  of  the  LOX  post  diameter.  Instability  mechanisms  in  coaxial  injectors  were  investigated 
experimentally  by  Glogowski  et  al.  [15,  16]  under  noncombusting  conditions.  Their  experiment  results  showed 
that  for  the  coaxial  injector  with  the  liquid  oxygen  (LOX)  post  recessed  into  the  fuel  annulus,  the  injector 
transitioned  into  a  condition  of  resonance  characterized  by  a  whistling  noise  and  significant  modification  to 
the  overall  structure  of  the  spray  due  to  the  strong  acoustic  coupling  between  injector  hydrodynamics  and 
spray  formation. 

Bazarov  [17,  18,  19]  studied  the  self-oscillation  phenomena  along  with  the  self-pulsation  mode  of  jet 
instability  in  coaxial  injector.  He  postulated  that  the  self-oscillation  occurs  when  the  gas-liquid  interaction 
forms  a  cavity  inside  nozzle,  leading  to  jet  swirling  around  the  nozzle  exit.  The  self-pulsation  of  the  liquid 
jet  mixed  with  the  gas  flow  depends  on  the  pressure  drop  at  liquid  and  gas  phase,  correlating  with  the  time 
of  liquid  propagation  through  the  injector  nozzle  [19]. 

In  this  report,  the  hydrodynamic  instability  of  the  coaxial  jet  in  the  recessed  region  is  presented  using  the 
three  dimensional  direct  numerical  simulation.  The  results  indicate  that  the  recess  length  to  injector  orifice 
diameter  has  a  significant  effect  on  spray  structure  over  the  velocity  and  gas  ratio  changes.  Unfortunately, 
none  of  experimental  data  for  exact  comparison  exists  due  to  very  small  spartial  and  temporal  scales  involved 
in  the  area.  However,  some  of  recent  experimental  data  provide  good  measure  in  a  macroscopic  point  of 
view  for  the  instability  frequency  as  explained  later  in  this  section. 

A  fully  unsteady,  three  dimensional,  two-phase  simulation  has  been  developed  utilizing  a  finite  volume 
implementation  of  the  Marker  and  Cell  discretization  method.  The  current  model  is  based  on  locally  ho- 
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Figure  2:  Sauter  mean  diameter  comparison  for  various  jet  speed  U  in  [m/s]. 


mogeneous  flow  (LHF)  assumption  in  which  the  relative  velocity  and  temperature  between  two-phases  are 
small  enough  in  comparison  to  variation  of  the  overall  flow  field  that  is  predicted.  Additional  constituitive 
relation  for  desity  field  has  been  implemented  in  order  to  provide  a  mechanism  of  solving  two-phase  flow 
with  a  single  phase  Navier-Stokes  equations  set.  This  fictious  “pseudo”  density  varies  in  amplitude  between 
the  liquid  and  gas  extream.  The  LHF  assumption  and  the  pseudo-density  implementation  allow  the  current 
model  to  handle  the  two-phase  flow  field  with  one  governing  equations  set  rather  than  to  compute  separate 
governing  equation  sets  for  each  flow  phase,  liquid  and  gas  in  this  case. 

The  code  runs  on  a  state-of-art  Beowulf  Linux  cluster  that  is  equipped  with  104  processors  and  fast 
ethernet  network.  One  run  usually  requires  12  to  24  processors  depending  on  mesh  discretization.  Since 
the  cluster  is  dedicated  to  the  modeling,  each  run  makes  use  of  nearly  100%  of  CPU  power  and  network 
bandwidth.  Even  with  this  superb  environment,  one  run  up  to  150,000  time  steps  takes  about  three  weeks. 

Parallel  processing  using  MPI  (Message  Passing  Interface)  has  been  implemented  in  order  to  run  the 
3-D  model  in  a  timely  manner.  For  the  present  3-D  atomization  modeling,  the  computational  domain  is 
split  up  in  axial  direction  for  the  desired  number  of  processors,  n.  While  each  processor  solves  a  flow  field 
of  subdivided  domain,  the  boundary  conditions  of  each  subdomain  are  transfered  to  neighboring  domain 
through  message  passing.  All  the  processors  involved  in  this  procedure  perform  the  same  calculation  in 
each  step  by  copying  the  original  task  to  other  processors  but  solving  different  subdomains  and  boundary 
conditions. 

The  Space  Shuttle  Main  Engine’s  (SSME)  coaxial  injector  has  been  chosen  for  baseline  modeling  case. 
The  injector  of  the  SSME  main  combustion  chamber  (MCC)  uses  liquid  oxygen  (LOX)  and  gaseous  hydrogen 
{OH 2)  as  an  oxidizer  and  fuel  respectively  in  fuel- rich  injection  environment  entering  the  main  combustion 
chamber. 

The  results  exhibit  characteristics  termed  ”  self-pulsation”  and  ” self-oscillation”  as  theorized  by  Bazarov  [17, 
18,  19].  Pulsations  are  evident  in  changes  in  the  flowrate  with  time,  while  oscillations  are  evidenced  by  az¬ 
imuthal  motion  of  the  central  liquid  core  about  the  annulus.  Prior  2-D  simulations  [20,  21]  have  shown 
similar  behavior,  although  the  azimuthal  motion  cannot  be  resolved  with  decreased  degree  of  freedom. 

The  overall  density  field  behavior  of  the  liquid  jet  for  conditions  roughly  equivalent  to  the  SSME  MCC 
injector  is  illustrated  in  Figure  4.  The  left  column  in  the  figure  depicts  density  contours  in  a  cutaway  view 
interpretation  of  the  motion  of  the  jet.  The  right  column  shows  the  density  contours  at  the  exit  plane. 
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Figure  3:  Effect  of  jet  speed  on  jet  surface  structure  of  Hoyt-Taylor’s  jet.  Due  to  larger  circulation  contained 
in  the  velocity  profile,  more  unstable  jet  structure  appears  in  the  lower  speed  case. 


The  resultant  highly  nonlinear,  quasi-periodic  oscillation  occurs  naturally  as  a  result  of  the  Kelvin- 
Helmholtz  instability  mechanism.  There  is  no  oscillation  in  the  inflow,  yet  strong  oscillations/pulsations 
develop  as  a  result  of  the  large  relative  dynamic  pressure  between  the  gas  and  liquid  streams.  The  develop¬ 
ment  of  self-pulsation  mode  of  the  liquid  jet  is  apparent  in  the  figure. 

One  interesting  bulk  measure  of  orifice  performance  is  the  massflow  delivered  at  the  orifice  exit  plane  as 
a  function  of  time.  This  quantity  can  be  computed  by  quadrature  of  the  density-axial  velocity  product  over 
the  exit  plane  area,  eigure  5  shows  the  time  history  of  mass  flow  for  the  density  field  computed  in  Fig.  4 
and  the  attendant  velocity  field.  The  unsteadiness  (self-pulsations)  of  the  jet  are  quite  evident  in  this  figure; 
massflow  variations  of  39%  about  the  mean  flow  are  apparent. 

The  spectral  content  of  the  mass  fluctuation  signal  provides  information  regarding  the  frequency  spectrum 
present.  Using  a  Fast  Fourier  Transform,  the  spectral  content  is  shown  in  Figure  6.  There  is  a  significant 
energy  content  over  a  range  of  dimensionless  frequencies  with  a  maximum  at  /  =  2.2074,  which  represent 
around  73,000  Hz  in  dimensional  units.  This  frequency  is  about  triple  the  frequency  at  which  liquid  is 
replaced  in  the  submergence  region. 

In  a  recent  study  of  Smith  and  Mayer  [22],  the  liquid  core  jet  oscillation  in  a  coaxial  injector  was  observed, 
and  the  frequency  of  flame  flickering  in  a  combustion  chamber  was  reported  as  near  7000  Hz,  which  falls  within 
the  same  frequency  band  of  the  SSME  main  burner  coaxial  injector’s  case.  The  same  velocity  and  density 
ratios  were  used  with  similar  geometry.  Even  though  their  coaxial  injector  was  not  exactly  the  same  one 
used  in  this  study,  their  experimental  study  provides  some  evidence  that  the  liquid  jet  oscillation/pulsation 
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plays  a  significant  role  in  combustion  response  within  the  chamber.  This  is  an  important  finding  requiring 
further  experimental  verification. 

The  three-dimensional  simulation  has  shown  a  good  predictive  capability  in  terms  of  coaxial  jet  behavior. 
This  allows  conducting  a  series  of  parametric  study  before  actual  injector  design  and  testing. 

3.3  Simulating  Droplet/ Wall  and  Droplet/Film  Collisions  Using  a  Level  Set 
Method 

A  numerical  method,  a  coupled  level  set  and  Marker  and  Cell,  is  developed  for  computing  axisymmetric, 
incompressible,  and  immiscible  two-phase  flows.  In  stead  of  using  marker  particles  to  track  the  free  surface, 
the  level  set  function  is  employed  to  “capture”  the  complex  interfacial  structure.  An  iterative  process  is 
devised  in  order  to  maintain  the  level  set  function  as  the  signed  distance  from  the  interface. 

As  a  base  line  case,  the  two-phase  fluid  code  is  implemented  for  simulating  zero  gravity  capillary  os¬ 
cillations  of  liquid  droplet  in  order  to  validate  its  capability  to  handle  the  surface  tension  effect.  Initially 
deformed  shapes  of  the  second  spherical  harmonics  are  to  oscillate  and  the  simulation  results  are  compared 
with  linearized  analytic  solutions  as  well  as  other  numerical  calculations.  We  investigated  the  viscous  effect 
on  oscillatory  damping  and  examined  the  role  of  surface  tension  on  oscillatory  period.  The  error  analyses 
relative  to  the  linearized  analytic  solutions  were  performed  by  the  successive  grid  refinement.  The  code’s 
numerical  accuracy  was  found  to  be  within  a  few  percent. 

An  impact  of  a  liquid  droplet  onto  solid  surface  and/or  shallow  liquid  layer  results  in  three  regimes: 
bouncing,  spreading,  and  splashing  depending  on  the  relative  importance  of  inertial,  viscous,  and  surface 
tension  effects.  Three  good  examples  among  various  experimental  works  on  a  single  droplet  impact  are 
chosen  to  be  simulated  numerically.  The  direct  comparisons  are  carried  out  between  the  experimental  and 
numerical  results.  Navier  slip  condition  is  to  be  assigned  as  a  velocity  boundary  condition  at  the  contact 
surface.  It  is  found  that  there  exists  a  criterion  to  distinguish  between  the  splashing  and  the  deposition 
events  in  terms  of  a  single  impact  parameter.  A  simple  dimensional  analysis  on  the  fluid  motion  during  the 
impact  and  the  deformation  process  is  performed  and  discussed. 


4  Professional  Activities 

The  efforts  outlined  in  the  previous  section  of  this  report  were  made  possible  by  two  grants  from  AFOSR.  A 
single  student,  Mr.  James  H.  Hilbing,  was  supported  under  the  base  grant  (F49620-94-1-0151).  In  addition, 
an  AASERT  grant  (F49620-93-1-0363)  was  utilized  to  support  Chris  A.  Spangler,  Mark  W.  Rutz,  Michael 
P.  Moses,  Ian  F.  Murray,  and  Kurt  Rump  (all  U.S.  citizens).  The  following  theses  were  written  as  a  result 
of  these  two  grants: 

Ph.P.  Dissertations 


Yoon,  S.  S.,  “A  Fully  Nonlinear  Model  for  Atomization  of  High-Speed-Jets”,  December,  2002. 

Kim,  B.,  “Study  of  Hydrodynamic  Instability  of  Shear  Coaxial  Injector  Flow  in  a  Recessed  Region”, 
December,  2002. 


A  list  of  journal  publications  (and  submissions)  associated  with  these  efforts  are  provided  in  the  following 
list.  Highlighted  items  (*)  have  been  attached  in  the  Appendices  of  this  report. 

Refereed  Journal  Publications  and  Submissions 


•  Yoon*,  S.  S.  and  Heister,  S.  D.,  “A  Fully  Nonlinear  Model  for  Atomization  of  High-Speed  Jets”, 
Engineering  Analysis  with  Boundary  Element ,  2002,  to  appear. 

•  Yoon,  S.  S.  and  Heister,  S.  D.,  “Analytic  Solutions  for  Computing  Velocities  Induced  from  Potential 
Vortex  Ring”,  AIAA  Journal ,  2002,  submitted. 
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•  Yoon*,  S.  S.  and  Heister,  S.  D.,  “Categorizing  Linear  Theories  for  Atomizing  Jets”,  Atomization  and 
Sprays ,  2002,  submitted. 

•  Yoon*,  S.  S.  and  Heister,  S.  D.,  “Modeling  Atomizing  Jet  due  to  Boundary  Layer  Instabilities”,  Physics 
of  Fluids ,  2002,  submitted. 

•  Kim*,  B  and  Heister,  S.  D.,  “Two-phase  Modeling  of  Hydrodynamic  Instabilities  in  Coaxial  Injectors” 
Journal  of  Propulsion  and  Power ,  2002,  submitted. 

•  Kim*,  B.  and  Heister,  S.  D.,  “Three-dimensional  Simulations  of  Flow  within  the  Recessed  Region  in  a 
Coaxial  Injector”,  In  Review,  J.  of  Propulsion  and  Power 

A  list  of  the  conference  papers  presented  in  association  with  work  under  these  grants  is  provided  in  the 
list  below. 

Conference  Papers  and  Presentations 


•  Yoon,  S.  S.  and  Heister,  S.  D.,  “A  Fully  Nonlinear  Primary  Atomization”,  15th  Annual  Conference  on 
Liquid  Atomization  and  Spray  Systems ,  pp.  36-40,  2002,  held  in  Madison,  Wisconsin. 

•  Yoon,  S.  S.  and  Heister,  S.  D.,  “A  Fully  Nonlinear  Primary  Atomization”,  38th  AIAA/ASME/SAE/ASEE 
Joint  Propulsion  Conference  and  Exhibit ,  AIAA  2002-4179,  held  in  Indianapolis,  Indiana. 

•  Kim,  B  and  Heister,  S.  D.,  “Two- phase  Modeling  and  Hydrodynamic  Instability  Study  of  Shear  Coaxial 
Injector  Flow”,  38th  AIAA/ASME/SAE/ASEE  Joint  Propulsion  Conference  and  Exhibit ,  AIAA  2002- 
3696,  held  in  Indianapolis,  Indiana. 

4.1  Technology  Transfer/Coupling  Activities 

Numerous  technology  transfers  have  occurred  during  the  period  associated  with  these  grants.  At  present, 
Professor  Heister  is  consulting  with  AFRL  officials  (Dr.  Tom  Hawkins  and  Dr.  Ron  Spores)  on  alternate 
ignition  schemes  for  high-performance  monopropellants  which  have  been  developed  at  AFRL.  Our  labora¬ 
tories  have  served  as  testing  facilities  for  Aerojet,  TRW,  Lockheed  Martin,  Pratt  &;  Whitney,  and  Boeing 
thereby  producing  useful  interactions  with  engineers  from  each  of  these  firms.  We  have  also  designed  en¬ 
gines/injectors  for  NASA  Dryden  Flight  Research  Center  and  a  small  company  (KB  Sciences)  which  was 
run  by  the  late  Dr.  Ron  Humble. 
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Figure  7:  Droplet  impact  deposition:  Droplet  impingement  with  low  impact  momentum  experiences  either 
(a)rebound  or  (b)spreading  depending  on  the  combination  of  Reynolds  and  Weber  numbers.  The  impact 
conditions  are  Re=100,  200x200  grids  in  dimensions  4x4  (a)  We=l,A£=.0005;  (b)  We=100,At=.0001. 
Gas/liquid  ratios  are  p'  =  10-3;//  =  10-2.  Gravitational  effect  is  ignored. 


Figure  8:  Limits  for  deposition/splashing  of  droplet 
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Figure  9:  Capillary  oscillation  of  liquid  droplet:  Released  from  a  static  deformation  of  second  (oblate)- 
spherical  harmonic  with  77  =  0.5,  a  zero  gravity  liquid  drop  is  oscillating  due  to  the  surface  tension.  This 
picture  spans  one  period  of  oscillation.  Re  —  100;  We  —  1;  200  x  100  grids;  pg/pt  =  Pg!  Pt  —  10 
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5  Appendix  A  -  Boundary  Element  Modeling 


Yoon*,  S.  S.  and  Heister,  S.  D.,  “A  Fully  Nonlinear  Model  for  Atomization  of  High-Speed 
Jets”. 
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ABSTRACT 

A  nonlinear  model  has  been  developed  to  assess  the  time  evolution  of  an  axisymmetric  liquid  jet  using  a 
boundary-element  method.  Vorticity  transported  from  the  boundary  layer  in  the  orifice  passage  to  the  free 
surface  is  modeled  using  a  potential  ring  vortex  placed  at  the  orifice  exit  plane.  The  vortex  strength  is 
uniquely  determined  from  the  Kutta  condition  and  information  regarding  the  boundary  layer  thickness  at 
the  orifice  exit  plane.  It  is  shown  that  primary  breakup  can  occur  even  in  the  absence  of  the  gas  phase. 
Using  a  secondary  stability  analysis  after  Ponstein  [1],  the  size  of  the  droplets  is  estimated  based  on  the  size 
of  the  ring-type  structures  shed  from  the  periphery  of  the  jet.  Computed  droplet  sizes  are  in  reasonable 
agreement  with  experimental  data,  although  turbulence  effects  obscure  some  comparisons. 

INTRODUCTION 

The  atomization  of  a  liquid  jet  is  one  of  the  most  fundamental  problems  of  two-phase  flow  and  has  received 
much  attention  due  to  the  large  number  of  practical  applications.  Since  the  formation  of  droplets  is  ultimately 
dictated  by  a  balance  of  capillary  and  inertial  forces,  numerical  methods  that  provide  high  resolution  of  these 
forces  perform  best  in  simulating  these  flows.  The  capillary  force  depends  on  local  surface  curvature,  which 
is  a  function  of  surface  shape  with  dependence  on  second  derivatives  of  local  surface  coordinates.  Resolving 
the  curvature  accurately  is  of  paramount  importance  in  these  problems.  For  these  reasons,  the  Boundary- 
Element  Method  (BEM)  is  uniquely  suited  for  atomization  modeling  in  that  the  optimal  placement  of  nodes 
on  the  gas/liquid  interface  provides  a  mean  for  maximizing  accuracy  of  surface  curvature  calculations. 

BEM  techniques  have  been  applied  to  a  wide  variety  of  free  surface/atomization  problems  including 
liquid  jets  [2,  3,  4,  5,  6,  7,  8,  9,  10],  droplets  [11,  12]  and  electrostatic  atomization  [13,  14].  Along  with  the 
difficulty  in  computing  surface  curvature,  the  inherent  nonlinearity  of  the  free-surface  condition  has  been 
addressed  by  many  researchers.  Several  proven  techniques  available  at  the  present  time.  In  the  liquid  jets 
field,  high-resolution  and  high  fidelity  BEM  solutions  for  low-speed  flows  are  now  available  and  applicable 
to  problems  in  chemical  engineering  and  inkjet  printing. 

Despite  these  advances,  the  atomization  process  increases  in  complexity  with  increasing  jet  speed.  For 
this  reason,  the  modeling  of  high-speed  jets  which  produce  small  droplets  of  interest  in  many  application, 
are  still  an  area  requiring  significant  efforts.  Prior  research  [3,  15,  16]  has  indicated  the  importance  of  the 
boundary  layer  structure  at  the  orifice  exit  plane  in  mapping  instabilities  just  downstream  of  the  injection 
point,  as  shown  in  Fig.  1.  The  present  study  focuses  on  this  issue  via  an  axisymmetric  simulation  that 
properly  accounts  for  the  presence  of,  and  vorticity  within,  this  boundary  layer.  The  following  section 
provides  a  description  of  the  model,  followed  by  convergence  studies  and  comparisons  with  experimental 
data. 


MODELING 

The  model  is  based  on  unsteady  axisymmetric  potential  flow  of  a  liquid  exiting  a  round  orifice  in  the  absence 
of  a  gas-phase  medium.  A  ring  vortex  is  employed  to  simulate  viscous  effects  associated  with  vorticity  in 
the  boundary  layer  formed  in  the  orifice  passage.  Carefully  controlled  experiments  have  shown  a  nearly 
axisymmetric  structure  during  the  early  stages  of  the  free-surface  instability.  Fig.  2  provides  a  schematic 
representation  of  the  geometry  and  nomenclature.  The  size  of  the  Rankine  vortex  [17]  is  defined  as  Rc,  as 
will  be  discussed  in  detail  in  a  later  section.  A  vortex  ring  of  strength  Tv  and  overall  radius  f  is  assumed  to 
lie  at  the  orifice  exit  plane.  A  computational  domain  is  represented  by  a  simple  cylindrical  column  of  length 
zL)  and  a  hemispherical  cap  is  selected  to  initialize  the  calculation.  Constant  nodal  spacing,  As,  is  employed 
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along  this  domain,  and  nodes  are  added  as  the  jet  issues  forward  from  the  orifice.  Fig.  3  compares  the  actual 
flow  condition  with  the  current  superposition  modeling.  Because  the  concentrated  vortices  at  the  filament 
vortex-ring  are  transported  to  the  free-surface,  the  jet  surface  is  unstable.  We  choose  the  liquid  density,  p, 
jet  average  exit  velocity,  U,  and  orifice  radius,  a  as  characteristic  dimensions  in  the  problem. 

The  formulation  of  the  BEM  starts  with  the  integral  representation  of  the  solution  of  Laplace’s  equation, 
V2^  =  0,  with  (j)  being  the  velocity  potential.  Following  Liggett  and  Liu  [18],  the  integral  form  for  this 
relation  is  given  by: 

ot<t>(ri)+  / 

Jn 

where  0(fl)  is  the  value  of  the  potential  at  a  point  f*,  ft  is  the  boundary  of  the  3D  domain,  and  G  is  the 
free  space  Green’s  function  of  Laplace’s  equation.  If  dQ,  is  a  3D  surface  eleilient,  ds  is  a  2D  surface  with 
dQ,  =  rdOds.  Then  the  governing  Eqn.  (1)  can  be  written  as: 


nr 


dfl  =  0 


(1) 


OL(j)(Ti) J  [(f>Dkern  Q^kern]  ds  —  0 


where  the  source  and  doublet  kernels  are: 
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and  the  auxilliary  quantities,  a,c,d  are  defined: 

c  =  (r  -  ri)(z  -  Zi)2 
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where  m  =  4r?',/a.  Here,  K(m)  and  E(m)  are  the  elliptical  integral  [19]  of  the  first  and  second  kind 
respectively,  and  nz>r  is  the  component  of  the  normal  vector  in  z  and  r  direction.  We  utilize  linear  elements 
for  the  velocity  potential, 

(5) 

\Sj+i  Sj  J  \Sj+ 1  sj  J 

with  an  analogous  form  for  the  normal  velocity,  d<f>/dn  =  q.  Substitution  of  Eqn.  (5)  into  Eqn.  (2)  yields 
the  following  Nnode  x  Nnode  matrix: 

[cd  +  D][4>]  =  [S\[q]  (6) 

where  a  arises  from  singular  contributions  to  the  doublet  terms  as  the  integral  passes  through  the  base  point 

Zj,Tj. 

Because  the  matrix  system  can  be  quite  large  (i.e.  Nnode  >  1000),  the  matrix  inversion  step  is  parallelized 
to  incorporate  multiple  processors.  Since  a  4th-order  Runge-Kutta  time  integration  scheme  is  employed,  the 
calculation  requires  solving  the  large  matrix  system  four  times  at  every  time  step.  In  addition,  the  size  of 
the  computational  domain  increases  in  time,  which  yields  a  larger  matrix  system  at  every  time  step  [20,  21]. 
The  Scalable  Linear  Algebra  Package  [22]  (SCALAPACK)  is  implemented  into  the  current  code  as  a  parallel 
solver.  Under  this  parallelization,  the  code  can  utilize  up  to  Np  number  of  processor.  We  have  used  Np  =  6 
because  little  variation  is  observed  for  the  computational  efficiency  after  Np  =  6  with  the  current  104  node 
LINUX-Cluster  networking  system.  Here,  a  tradeoff  exists  between  individual  node  speed  (we  use  1.2  GHz 
Athalon  chips)  and  nodal  communication  time.  Because  the  nodes  are  quite  fast  in  performing  the  required 
calculations,  communication  time  tends  to  dominate  for  Np  >  6. 

A  typical  [A][X]  =  [B]  non-parallel  solver  based  on  LU-Decomposition  is  available  in  Numerical  Recipes  [23]. 
When  the  matrix  [A]  and  vector  [B]  are  given,  the  matrix  [X]  can  be  obtained  using  the  LU-Decomposition. 
The  parallel  [A][X]  =  [B]  solver  SCALAPACK  is  also  based  on  the  LU-Decomposition  algorithm.  However, 
one  additional  step  is  required  for  SCALAPACK:  the  partition  of  [A],  [X],  and  [B]  into  the  Np  process  grid. 
Unfortunately,  SCALAPACK  does  not  perform  the  matrix  partition  autonomously.  Thus,  we  have  developed 
the  algorithm  for  the  grid-partition.  Considering  Np  —  6,  the  physical  matrix  [A]  is  transformed  into  the 
2x3  process  grid,  as  shown  in  Fig.  4.  Parameters  that  define  this  transformation  are  MB  x  NB  for  the  size 
of  the  submatrix,  and  NPROW  x  NPCOL  for  the  number  of  processors  (see  Fig.  4).  The  particular  parallel 
solver  used  for  our  calculation,  called  ‘pdgesv’  [22],  can  handle  the  fullly  dense  matrix,  [A].  If  the  matrix 
[A]  is  tridiagonal,  and  ‘pddtsv’  is  used  for  the  efficient  calculation.  In  addition,  ‘pddbsv’  solves  a  general 
band  system  [A][X]  =  [6].  It  should  be  noted  that  the  answer  [X]  is  splitt  into  crow’  processors.  As  shown  in 
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Fig.  5,  we  have  combined  the  answer  [X]  in  the  (0,0)  processor  and  broadcast  [X]  to  other  processors  using 
the  MPI  routine. 

We  have  constructed  the  sample  calculation;  Laplace  equation  with  the  Dirichlet  boundary  conditions  at 
the  side  walls  and  thus  [A]  is  a  banded  system  and  ‘pdgesv’  is  used.  The  sample  code  can  be  downloaded  at 

http://roger.ecn. purdue. edu/~  yoons/SCAL.tar.gz 

Flow  physics  at  the  smaller  length  scale  can  be  captured  using  a  finer  grid  under  the  parallelization.  The 
length  scale  of  our  interest  is  scaled  by  the  momentum  thickness,  5 2.  A  complete  calculation  takes  about  one 
to  two  weeks:  this  is  long  enough  real  time  to  obtain  all  the  information  about  the  droplet  characteristics  of 
our  particular  interest. 


Superposition 

Contributions  from  the  vortex  ring  can  be  obtained  through  the  principle  of  superposition  for  potential  flow. 
Since  Laplace’s  governing  equation  is  linear,  we  may  superpose  the  bulk  potential  flow  with  the  flow  due  to 
the  vortex  ring  (i.e.  <j>v  or  tyv), 

</>£  =  4>  +  <t>v  ®  +  Vv  (7) 

where  <j)  and  Psi  represent  the  velocity  poential  and  the  streamfunction,  respectively.  Here  (  )£  denotes 
the  total  solution  of  the  jet  flow.  The  flow  of  the  vortex  ring  can  be  obtained  either  by  applying  the 
Cauchy- Riemann  condition  to  the  stream  function  [24],  or  by  direct  evaluation  using  Biot-Savart  law  [25]. 
Traditionally  [24,  26,  17],  the  velocities  induced  from  this  flow  are  computed  from  the  stream  function, 


$(z,r)  = 


K{m)  -  —r =E{m ) 

■y  771 


(8) 


where  m  =  m(r,z).  Here,  I\,  is  the  strength  of  the  vortex  and  K (m)  and  E(rn)  are  the  complete  elliptic 
integrals  of  the  first  and  second  kind,  respectively.  The  velocities  can  be  determined  in  a  standard  fashion 
by  differentiating  'f, 


uv  =  - 


13$ 


r  dr 


13$ 

Vv  = - -X- 

r  az 


(9) 


Alternatively,  the  Biot-Savart  law  [27,  28,  25]  can  be  evaluated  directly  for  the  filament  vortex- ring  solution, 


uv 


x  (f  -  fj) 
|  f-fi  |3 


ddi 


(10) 


where  dl  is  the  path  enclosed  by  the  vortex  ring.  The  entire  explicit  solution  of  the  filament  vortex-ring  has 
recently  been  developed  [29]. 


Free  Surface  Boundary  Conditions 


The  unsteady  Bernoulli  equation  provides  the  boundary  condition  along  a  free  surface  [30].  This  condition 
provides  a  connection  between  the  inertial,  hydrostatic,  and  capillary  forces  at  the  interface.  Because  the 
surface  curvature,  k,  depends  nonlinearly  on  the  surface  shape,  the  overall  expression  is  nonlinear.  Using  the 
nondimensionalization  described  previously,  the  appropriate  dimensionless  form  of  the  Bernoulli  equation  is 
given  by: 


^P’  +  Wl-Wi 


z  —  0 


(ii) 


where  Pg  is  the  dimensionless  gas-phase  pressure  (assumed  to  be  zero  in  the  present  studies),  and  We  and 
Bo  are  the  Weber  and  Bond  numbers, 


We  = 


pU2a 

a 


(12) 


where  a  is  the  surface  tension  and  g  is  the  gravitational  constant.  As  shown  in  Fig.  2,  the  vortex  ring  is 
located  at  the  nozzle  exit  in  order  to  model  the  vorticity  present  in  the  actual  flow.  While  the  singularity 
arises  at  the  vortex  ring,  the  solution  at  the  vortex-ring  is  also  singular  and  that  of  the  computational 
nodes  near  the  vortex-ring  is  not  physical.  In  an  effort  to  overcome  the  problem,  we  have  separated  the 
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governing  equation  for  potential  flow  from  the  total  solution;  the  time  derivative  of  the  vortex-ring  is  zero 
(i.e.  d<t>v/dt  =  0)  since  its  strength  (or  circulation)  is  assumed  constant  and  its  location  is  fixed.  Thus  the 
time  derivative  of  the  Eqn.  (7)  is  given  by: 

d(/>t  _  (ty 
dt  dt 

Thus  Eqn.  (11)  becomes: 

dfi  1,„  ,  ,o  k  Bo 

~dt  2  +Pg+We~W^Z~ 

The  transformation  from  the  Eulerian  (i.e.  d/dt)  to  the  Lagrangian  (i.e.  D/Dt)  is  required  as  the  compu¬ 
tational  nodes  move  respect  to  time,  t: 


(13) 

(14) 


Combining  Eqns.  (14)  and  (15)  yields: 


2£  =  v*.v*-i|  +  ^ 

Substitution  of  Eqn.  (7)  into  Eqn.  (16)  (for  Wfi  term)  yields: 

^  =  Vfit  •  V(</>t  -  fiv)  -  ^|V(£t|2  -Pg--^-  + 

1,„,  l9  ^  k  Bo 

=  •  Vfit  -  V<f>t  ■  Vfiv  -  -Pg-  +  yy-Z 


(15) 


(16) 


(17) 


1 


Bo 


=  -|V^|2-V^-V^-Pg-  — +  ^ 

Since  differential  operator  (V)  is  a  linear  function,  superposition  theory  holds  for  the  velocity  as  well: 

ut  =  u  +  uv  Vt=v  +  vv  (18) 


Combining  Eqns.  (17)  and  (18),  we  obtain  more  explicit  Bernoulli  relation  for  fi: 


P0  _  1  |„,  |2 


Bo 


—  —  -|ut|  Ut  ■  Uv  Pg  + 


Dt 


We 


(19) 


Now  that  we  have  separated  the  fi  from  the  general  solution  fit,  the  entire  time  marching  of  the  numerical 
algorithm  is  related  to  fi  only.  The  effect  of  the  vortex  ring  comes  into  play  only  in  Eqn.  (18).  Bernoulli’s 
Eqn.  (19)  is  the  ‘nonlinear’  boundary  condition  for  the  Laplace  equation  and  it  is  marched  in  time  using  a 
4th  order  Runge-Kutta  time  integration. 

The  curvature  of  the  highly-distorted  surface  is  determined  with  4th-order  accuracy  [31,  32].  We  have 
used  the  curvature  definition  from  Smirnov  [33]: 


The  derivatives  are  evaluated  using  centeral  differencing  with  the  exception  of  the  ends  of  the  domain  where 
forward  or  backward  derivatives  are  applied  [32]. 

The  location  of  nodes  on  the  free  surface  (i.e.  z  and  r)  is  advanced  by  integrating  the  respective  velocity 


components  in  time: 

Dz  d(j> 

Dr 

_  d(f> 

(21) 

Dt  dz 

~Dt 

dr 

where 

dfi 
dz  = 

96  .  ^ 

=  —  cos  p  -  q  sin  p 

9s 

dfi 

dr 

96 

—  —  sin  p  4-  q  cos  p 

9s 

(22) 

where  q  =  dfi/dn  is  the  normal  derivative  of  phi.  The  surface  slope,  /?,  is  assumed  to  be  given  by  the  slope 
of  the  parabola  at  the  middle  node  [31]. 
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Smoothing 

Nodes  are  repositioned  along  the  distorted  surface  using  cubic  splines  [34],  and  nodes  can  be  added  (due  to 
fluid  exiting  the  nozzle)  or  removed  (due  to  atomization  events)  without  user  intervention.  The  current  ‘high- 
speed’  atomization  simulation  is  more  susceptible  to  numerical  instability  than  the  low-speed  atomization 
(i.e.  varicose  breakup  simulation  with  smooth  surface)  due  to  complex  surface  shape.  In  addition,  the 
‘necking’  region  where  a  droplet  is  pinched-off  experiences  a  relatively  high  velocity  (i.e.  node  velocities 
reache  2  ~  7  times  that  of  the  jet  speed)  and  therefore  ‘node-crossing’  or  ‘not-simply-connected  domain’ 
sometimes  occurs.  Thus  numerical  smoothing  is  crucial  to  prevent  the  numerical  instability  [35]. 

Following  the  definition  for  filter  function  and  its  transfer  function  by  Spyropoulos  and  Blaisdell  [36]: 

7  _  „  /j-n  +  /j+n  ^  ( 2irkn\ 

ft  -  a* - o -  Gk=l^  a"cos  (  -jr  )  (23) 

n=0  n=0  \  iV  / 

where  /  is  the  filtered  function,  /  is  the  function  prior  to  the  filtering,  an  is  the  coefficient  listed  in  Table  1 
for  N  —  2  4,  N  is  the  number  of  points  used,  Gk  is  the  filter  transfer  function,  and  k  is  the  wavenumber. 

Table  1  provides  the  coefficients  for  the  various  explicit  filter  functions  developed  by  researchers  [37,  11, 
35,  38,  36].  In  Fig.  6,  the  transfer  filter  function,  G*,  is  plotted  as  a  function  of  the  wavenumber,  2k/ N. 
The  main  purpose  of  using  filter  is  to  suppress  the  relatively  large  wavenumbers  (or  small  wavelengths). 
The  filters  by  Hilbing  et  al  [37]  and  Lundgren  and  Mansour  are  not  suitable  because  they  do  not  suppress 
large  wavenumbers  sufficiently.  The  filter  by  Spyropoulos  and  Blaisdell  [36]  is  not  suitable  either  as  it 
amplifies  some  wavenumbers.  Since  the  3-pt  trapezoid  filter  damps  out  the  small  wavenumbers  more  than 
other  comparable  filters,  it  may  be  too  diffusive  for  the  current  use.  While  we  hope  to  damp  out  the  large 
wavenumbers  and  leave  the  small  wavenumbers  as  they  are,  we  have  chosen  the  filter  function  by  Longuet- 
Higgins  and  Cokelet  [35];  this  filter  is  designed  to  eliminate  the  odd-even  mode  in  the  function  which  contains 
the  highest  frequency,  known  as  the  Nyquist  Rate  [39]. 


Table  I:  Coefficients  for  the  various  explicit  filter  functions 


N 

a0 

ai 

a3 

3-pt  Trapezoid 

2 

0.5 

0.5 

— 

— 

5-pt  Hilbing  et  al.  [37] 

3 

0.9625 

0.05 

-0.0125 

— 

5-pt  Lundgren  and  Mansour  [11] 

mm 

0.94 

-0.02 

— 

5-pt  Longuet-Higgin  and  Cokelet  [35] 

mm 

0.625 

0.5 

-0.125 

— 

7-pt  Lele  [38] 

mm 

0.5 

— 

-0.0625 

7-pt  Spyropoulos  and  Blaisdell  [36] 

KM 

0.5 

0.6744132 

— 

Vorticity  Centroid 

The  centroid  of  vorticity  of  the  viscous  flow  is  regarded  as  the  center  of  the  vortex-ring.  The  definition  of 
the  centroid  of  the  vorticity,  weighted  in  the  radial  direction  is  given  by: 


r  = 


rujdr 
fr=0  udr 


(24) 


and  the  vorticity,  a;,  is  defined  as: 


dv  du 
dz  dr 


(25) 


where  u  and  v  are  the  velocities  in  the  axial  and  radial  direction,  respectively.  Assuming  dv/dz  ~  0,  i.e. 
parallel  flow  at  the  orifice  exit  plane,  the  centroid  can  be  written  as: 


f  = 


ru(r-l) 

Ju(r—0) 

ru(r=l) 
Ju(r— 0) 


rdu 

du 


(26) 
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Substituting  the  definition  of  the  displacement  thickness,  Jj,  [40]  into  Eqn.  (26)  and  applying  integration  by 
parts  gives  the  following  result: 

f  =  1  —  <5i  (27) 


where 


(28) 


can  be  approximated  using  a  Navier-Stokes  solution  of  the  internal  flow,  or  appropriate  analytical  methods 
such  as  Blasius  solution  [40]  for  a  flat-plate  or  the  Thwaites’  [41]  equation  for  a  converging/diverging  nozzle. 
The  vortex  strength  T„  is  defined  as  the  circulation  around  any  path  enclosing  the  vortex-ring,  namely: 


r„  = 


(29) 


where  u  is  the  internal  flow  velocity  of  the  injector  orifice  and  I  is  the  integration  path.  The  integration 
is  taken  in  a  counterclockwise  direction  around  the  path.  Assuming  du/dz  =  0  and  a  no-slip  boundary 
condition  at  the  wall,  the  surface  integral  in  Eqn.  (29)  can  be  written  as  follows  for  the  path  we  have  chosen: 

T„  =  Az  (30) 


where  A z  is  the  length  scale  which  is  comparable  to  the  most  dominant  wavelength,  A.  T„  is  always  positive 
and  will  induce  counterclockwise  motion  (this  is  based  on  the  upper  half  of  the  flow  going  from  left  to  right). 
Eqns.  (27)  and  (30)  uniquely  determine  the  location  and  strength  of  the  vortex  from  first  principles.  No 
additional  calibration  constants  are  used  in  the  formulation. 

It  is  possible  to  set  Az  =  A,  predicted  by  Brennen’s  result  [42],  since  the  circulation  causes  the  axisym- 
metrically  disturbed  wavelength  observed  at  the  nozzle  exit.  The  A  observed  in  Hoyt  and  Taylor’s  case  [16]  is 
a  function  of  the  momentum  thickness  52,  scaled  by  the  parameter  7  =  0.175.  Presuming  a  high  contraction 
ratio  of  the  nozzle  reduces  the  turbulence  fluctuation,  Hoyt  and  Taylor  [16]  assumed  a  laminarized  flow  over 
a  flat  plat  and  therefore  they  utilized  the  Blasius  [40]  solution  to  approximate  the  momentum  thickness, 


Another  possible  choice  for  the  circulation  approximation  is  to  use  the  Kutta-condition  [30]  from  the 
potential  flow  theory.  It  is  known  that  the  velocity  of  the  potential  flow  at  a  sharp  corner  is  infinite  and 
therefore  it  is  not  physically  possible.  One  remedial  treatment  is  that  a  stagnation  point  is  forced  at  the 
sharp  corner,  known  as  the  ‘Kutta-condition’.  While  the  gradient  of  the  inflow  velocity  of  the  injector  orifice 
is  zero  due  to  the  nature  of  potential  flow,  one  may  adjust  the  strength  of  vortex-ring  r„  to  be  such  that  the 
sum  of  potential  velocity  and  vortex  induced  velocity  in  the  axial  direction  to  be  zero  at  the  sharp  corner, 
this  linear  relation  was  applied  and  solved  to  find  the  r „ .  The  closer  the  vortex-ring  is  to  the  wall  (i.e.  the 
higher  f),  the  stronger  the  vortex  induced  velocity  will  be.  In  order  to  satisfy  the  stagnant  velocity  at  the 
sharp  corner  with  constant  potential  inflow  velocity  of  1.0,  r„  must  decrease  as  f  increases.  It  is  shown  in 
Fig.  7  that  the  strength  of  the  vortex-ring  condition  decreases  as  f  increases.  One  may  utilize  the  cubically 
interpolated  equation  of  the  Kutta-condition  in  Fig.  7 . 

The  third  possible  choice  is  to  adjust  the  strength  to  be  such  that  the  vorticity  integral  over  the  radial 
direction  of  the  real  flow  matches  that  of  the  vortex-ring: 


rr=  1  pr—  1 

/  udr  =  /  LJvdr  (32) 

J T= 0  Jt— 0 

It  is  known  that  uiv  is  infinite  at  the  ring  and  zero  everywhere  else.  However,  the  change  in  the  induced 
axial- velocity  over  the  radial  direction  is  finite  (i.e.  duv/dr).  Thus  we  again  assume  a  parallel  flow  for  the 
real  flow  (i.e.  dv/dz  =  0).  Then  Eqn.  (32)  becomes: 


ru(r=  1)  rUv[T=L) 

/  du=  duv 

Ju(r— 0)  Juv(r= 0) 


uv(r= 1) 


(33) 


This  is  essentially  matching  A uv  to  be  U  and  will  yield  the  positive  IV  The  quadratic  interpolation  of  this 
approximation  is  shown  in  Fig.  7. 
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The  variation  of  rv  as  a  function  of  the  centroid  of  the  vorticity  f  is  presented  in  Fig.  7.  It  is  seen  that 
the  rv  used  is  lower  than  that  of  the  Kutta-condition.  The  ratio  of  their  Tv  value  varies  from  2.8  to  1.0 
as  a  function  of  r.  It  is  observed  that  their  difference  in  Tv  decreases  as  f  increases  and  eventually  yields 
difference  of  approximately  3%  from  one  to  the  other  at  about  f  =  0.99.  The  reason  for  such  a  behavior  is 
that  the  model  accuracy  improves  when  the  boundary  layer  is  ‘thin’. 

In  summary,  any  choice  of  the  three  methods  is  applicable  for  the  circulation  approximation  not  only 
because  differs  little  at  high  f  but  also  because  Tv  by  itself  does  not  change  the  droplet  size  significantly. 
However,  one  needs  to  use  the  same  approximation  consistently  when  parametric  studies  are  performed.  We 
have  used  the  approximation  using  Brennen’s  theory  for  the  case  presented  in  the  result  section. 


Post  Processing  Formulation  for  Droplets 

The  addition  of  the  Biot-Savart  Law  to  the  inviscid  jet  of  the  BEM  formulation  is  expected  to  cause  instability 
at  the  free  surface  that  eventually  forms  a  series  of  toroidal  ligaments  pinched  off  from  the  main  body  of  the 
jet.  Here  the  formulations  for  the  cross-sectional  area,  centroids  of  the  area,  volume,  and  the  velocity  of  the 
droplets  are  developed. 

Using  Gauss’s  divergence  theorem,  we  have  transformed  the  surface  integral  to  a  line  or  contour  integral 
in  order  to  compute  the  cross-sectional  area  of  a  pinched-off  ring: 


A‘WcrdZ  +  IcZdT} 

which  was  discretized  using  the  trapezoid  Rule, 


(34) 


1  N 

A  ~  7  +  ri+i)(2i+i  -  ZJ )  -  (zi  +  zj+i)(rj+i  -  rj)\ 


3=1 


(35) 


where  N  is  the  last 
defined  as: 


(or  maximum)  node  that  closes  the  loop.  The  centroids 


= 


fs 

A 


of  the  cross-sectional  area  are 


(36) 


This  is  also  discretized  using  the  trapezoid  Rule.  Similar  to  the  previous  approach,  we  have  found  the 
cross-sectional  area  can  be  expressed  as: 


1 

2c  =  1^4  X!  [(rJ  +  r3+i)(zj+i  -  Zj)  -  ( Zj  +  zj+1)2(rj+1  -  Tj)]  (37) 

3= 1 
1  N 

r°  ~~  12 A  ~  zj)  ~~  izj  “b  2j+i)(r^+1  —  Tj) j  (38) 

3= i 

The  theorem  of  Pappus-Guldinus  [43]  relates  a  volume  of  revolution  to  its  generating  cross-sectional  area: 


V  =  2-kAtc  (39) 

For  the  velocities  of  a  droplet,  we  have  weighted  both  velocities  uD  (axial)and  vD  (radial)  in  the  radial 
direction.  The  formulation  is  as  follows: 


Up  = 


(40) 


Modeling  Secondary  Instability 

While  the  current  model  is  based  on  an  axisymmetric  formulation,  actual  primary  atomization  is  a  3- 
dimensional  phenomenon  (see  Fig.  8-(a)).  The  model  result  of  pinch-off,  as  shown  in  Fig.  8-(b),  is  not  a 
droplet  but  a  vortex-ring  with  a  significant  amount  of  circulation  around  the  ring  surface.  The  circulation 
around  the  ring  surface  is  large  enough  to  cause  instability  in  the  circumferential  direction.  In  reality,  this 
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is  the  secondary  instability  which  occurs  before  the  vortex-ring  pinch-off.  The  current  model  assumes  that 
droplets  are  formed  from  a  secondary  instability  on  annular  ligaments  shed  from  the  periphery  of  the  jet. 
This  amounts  to  a  decoupling  of  primary  and  secondary  instability  which  permits  the  axisymmetric  analysis 
of  the  jet  itself. 

Ponstein[l]  investigated  the  linear  stability  of  a  liquid  column  with  circulation  Tr  and  radius  aT.  Pon- 
stein’s  result  is  utilized  to  predict  the  dominant  wavelength  k  in  the  circumferential  direction: 


u 


2 


(41) 


where  Wer  -  pU2ar/cr.  Note  that  the  ring  radius,  ar,  is  the  non-dimensional  variable.  This  expression  is 
solved  to  determine  the  k  =  k^ax  value  attributed  to  the  maximum  growth  rate,  uj  for  a  given  ring  geometry 
and  circulation.  Since  Ponstein’s  analysis  was  conducted  for  a  liquid  column,  we  assume  that  the  thickness 
of  the  ring-shaped  ligaments  is  much  less  than  the  nozzle/jet  radius  (i.e.  ar  «  a).  This  assumption  is 
confirmed  from  ligament  sizes  produced  in  the  calculations.  Fig.  9  illustrates  how  Ponstein’s  equation  is 
applied  to  the  vortex-ring  (annular  ligament)  with  circulation  Tr. 

At  the  event  of  vortex-ring  pinch-off,  all  information  about  the  ring  is  collected  (i.e.  volume,  centroids, 
etc).  Thus  we  can  calculate  the  circumferential  length  of  the  ring,  l  =  2tt rc.  The  most  dominant  wavelength 
(A  =  2n /kmax)  which  corresponds  to  the  maximum  growth  rate  is  known  from  Ponstein’s  Eqn  (41).  Thus 
the  number  of  droplets  per  ring  can  be  estimated  (i.e.  N&  =  //A).  The  volume  of  the  ring  (i.e.  V)  is  known 
and  that  of  a  droplet  can  be  approximated  using  the  definition  of  a  sphere  volume,  Vd  =  V/ND  =  f  7r  (^)  . 
Thus  a  droplet  diameter  ( D )  can  be  estimated. 


RESULTS 

Grid  Convergence  Study 

Hoyt  and  Taylor’s  case  is  used  for  the  grid  convergence  check  (i.e.  We  =  19057,  r  =  0.99,  and  Tv  =  0.139). 
As  is  the  grid  spacing  for  the  BEM  node.  While  Hilbing  [32]  mentioned  that  As  =  0.300  is  fine  enough  to 
resolve  the  low  speed  “Rayleigh  breakup”  where  waves  are  of  length  comparable  to  the  orifice  diameter,  a 
much  finer  grid  resolution  is  required  for  high  speed  atomization  where  the  wavelengths  are  comparable  to 
the  boundary  layer  thickness  at  the  orifice  exit.  For  this  reason,  the  grid  resolution  for  the  present  studies 
taxes  the  current  computational  capabilities  of  even  advanced  Linux-based  compute  clusters.  In  Fig.  10, 
it  is  shown  that  the  axial  location  for  the  first  ring  pinch-off  is  reasonably  insensitive  to  mesh  spacing  for 
As  £  0.030.  However,  grid  function  convergence  studies  indicate  that  a  smaller  mesh  spacing  is  required  for 
the  accurate  prediction  of  the  droplet  characteristics  in  the  atomization  regime.  About  1000  ~  3000  droplets 
were  collected  for  each  run  for  statistically  reliable  data  and  results  for  drop  statistics  are  shown  in  Table  2. 
The  Sauter  Mean  Diameter,  SMD,  (drop  whose  diameter  replicates  the  average  surface  area  ofdrops  in  the 
population)  is  the  most  frequent  measure  used  in  the  atomization  field.  Table  2  shows  that  the  N u  per  ring, 
its  Standard  Deviation,  and  SMD  are  converged  to  a  reasonable  accuracy  at  As  =  0.012.  Thus  we  have 
used  As  =  0.016.  In  addition,  the  standard  deviation  of  ND/ring  does  not  change  much  after  As  ^  0.016. 
Similarly,  the  time-averaged  droplet  velocities  (i.e.  uo  and  vp)  do  not  change  much  either  for  this  mesh 
spacing.  This  result  also  validates  the  pinch  criteria  employed  for  a  ligament  breakup:  a  pinch-off  is  assumed 
when  the  distance  between  binary  nodes  is  less  than  a  certain  tolerance,  e.  The  range  of  the  tolerance  is 
20  ~  70%  of  the  mesh  spacing.  The  ligament  size  is  also  insensitive  to  the  range  of  the  pinch-off  criteria. 

It  is  uncertain  when  to  stop  the  simulation  since  the  jet  can  grow  indefinitely  depending  on  injection 
conditions.  For  the  simulations  conducted  to  date,  the  time  required  for  the  first  pinching  event  is  typically 
around  i  «  1.7.  We  found  that  collecting  about  300  to  400  rings  provides  stastically  reliable  data.  This 
would  give  roughly  1000  droplets.  Thus  we  typically  stop  our  calculation  at  about  t  ~  5.0.  Table  3  shows 
that  the  solution  is  insensitive  to  what  for  t  >  4.0. 

Effect  of  Rankine  Vortex  Size  and  Initial  Jet  Length 

In  Fig.  1,  vortices  induce  motion/instability  near  to  the  nozzle  exit  (i.e.  axisymmetrically  disturbed  waves) 
and  eventually  cause  the  jet  to  break  up  into  turbulent  flow.  However,  the  flow  at  the  nozzle  exit  is 
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Table  2:  Grid  convergence  test 


N  o/rinp 

Stnd.  Dev. 

SMD/rf 

MUM 

10.07 

4.68 

10.87 

8.71 

EZESI 

14.20 

11.33 

12.46 

8.46 

0.770 

■h»i« 

10.44 

6.65 

gjl^u 

10.87 

6.44 

ISEOiI&fl 

EXEM 

Table  3:  Effect  of  Calculation  Length  on  Drop  Statistics 


t 

SMD/cf 

No 

N  d /ring 

Stnd.  Dev. 

2.0 

0.0628 

88 

7.97 

3.329 

3.0 

0.0623 

856 

10.44 

6.665 

4.0 

0.0635 

2582 

11.08 

7.936 

5.0 

0.0655 

5132 

11.48 

7.985 

nearly  laminar  flow  as  the  flow  is  laminarized  through  highly  contracted  nozzle  geometry  under  a  favorable 
pressure  gradient  [16].  Since  the  filament  vortex-ring  is  located  exactly  at  the  nozzle  exit,  the  computational 
nodes  near  the  nozzle  exit  are  induced  with  the  greater  motion.  This  seems  contradictory  to  the  observed 
laminarized  flow  as  shown  in  Fig.  1.  In  reality,  it  takes  some  time  and  distance  for  the  rollup  motion  to 
develop  and  therefore  the  relaxation  length  is  present  regardless  of  the  flow  regime.  For  this  reason,  a  cutoff 
for  the  superimposition  method  of  the  filament  vortex-ring  is  introduced  using  stationary  Rankine  vortex 
model  [17].  The  size  of  Rankine  vortex  (i.e.  Rc  in  Fig.  2),  whose  center  is  located  at  the  upper  corner  of 
the  nozzle  exit,  has  little  effect  on  droplet  size  and  thus  this  is  the  parameter  that  can  be  set  at  the  users 
convenience.  As  shown  in  Fig.  11,  little  variation  in  the  axial  pinch-off  location  is  observed  for  Rc  <  0.4.  We 
have  set  Rc  =  0.3  so  that  the  computational  nodes  at  the  near  nozzle  exit  are  not  affected  by  the  induced 
motion.  This  is  essentially  setting  the  Rankine  vortex  size  to  be  the  relaxation  length,  Rc&lr.  It  should  be 
noted  that  the  relaxation  length  can  be  scaled  with  the  nozzle  length,  l  [44,  45]. 

The  simulation  begins  with  an  initially  assumed  jet  shape  as  shown  in  Fig.  2.  The  size  of  the  initial  shape 
seems  to  have  some  effect  on  the  first  pinch-off.  The  larger  the  shape,  the  earlier  the  motion/instability  is 
induced.  However,  the  initial  shape  size  seems  to  have  little  effect  on  the  axial  pinch-off  location  as  shown 
in  Fig.  12. 

Comparison  with  Experiment 

The  complete  simulation  of  the  Hoyt  and  Taylor’s  jet  [46]  is  shown  in  Fig.  13.  The  jet  structure  is  initially 
assumed  to  be  a  simple  cylinder  with  a  hemispherical  tip  as  shown  in  Fig.  2  and  its  evolution  is  simulated 
via  time  integration.  The  simulation  is  completed  at  t  =  5.0.  A  slight  ‘swelling’  is  observed  at  t  =  1.0  and 
the  fluctuation  of  the  jet  surface  is  seen  at  t  >  2.0.  The  velocities  induced  by  the  bound  vortex  are  large 
enough  to  penetrate  the  jet  surface  and  it  results  in  the  primary  atomization.  It  should  be  noted  that  most 
liquid  ligaments,  pinching  from  the  jet  surface,  are  in  the  ‘rollup’  motion  in  counterclockwise  direction  while 
the  mean  velocity  of  the  ligament  is  in  the  streamwise  direction.  Similar  structures  are  noted  in  Fig.  8- (a)  in 
a  closeup  view  of  the  Hoyt  and  Taylor  experiment.  The  counterclockwise  rollup  motion  is  a  strong  evidence 
that  the  boundary  layer  instability  is  the  fundamental  cause  of  the  primary  atomization.  The  mean  velocity 
of  most  droplets  are  in  streamwise  direction  as  the  droplets  motion  propagates  along  with  the  main  jet 
stream,  the  most  dominant  convective  source. 

It  was  mentioned, in  the  previous  section,  that  the  model  takes  advantage  of  Ponstein’s  Eqn.  (41)  to 
model  the  instability  of  a  pinch-off  vortex-ring.  Using  the  Eqn.  (41),  the  number  of  waves  or  droplets  (ND) 
per  ring  is  predicted  and  plotted  as  a  function  of  the  circulation,  Tr.  Their  relationship  is  parabolic  and 
the  least  square  fit  is  available  in  Fig.  14.  Thus  if  Tr  is  known,  a  rough  estimation  of  N&  can  be  produced. 
Fig.  15  illustrates  the  relation  between  the  number  of  circumferential  waves  and  the  core  size  of  a  pinch-off 
vortex- ring.  As  the  size  increases,  more  waves  appear.  This  is  exactly  the  opposite  phenomenon  as  for  the 
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elliptical  instability  of  Widnall  [47]  and  Sullivan:  they  observed  more  waves  when  the  core  size  decreases.  It 
should  be  noted  that  the  mechanism  of  the  surface  tension  driven  instability  (i.  e.  liquid  vortex-ring  in  air) 
is  different  from  the  shear  layer  driven  elliptical  instability  [48,  49,  50,  51],  (i.e.  liquid  into  liquid  or  gas  into 
gas). 

It  is  well  known  that  the  droplet  size  varies  significantly  within  the  atomization  regime.  Wu  et  al  [52] 
reported  the  droplet  size  variation  with  U  for  the  turbulent  water  jet  into  air.  Hoyt  and  Taylor  s  experiment 
had  been  carried  out  for  the  Bernoulli  pressure  A P  <  60  psi;  no  result  with  higher  A P  is  reported  [46, 16,  53]. 
However,  we  hypothesized  the  increase  in  A P  for  the  Hoyt  and  Taylor  jet.  The  result  is  taken  the  jet  speed 
up  to  17  =  40  m/s  which  corresponds  to  A P  ~  116  psi  or  slightly  higher  in  order  to  account  for  some  pressure 
loss  within  the  nozzle.  The  final  result  for  the  Hoyt  and  Taylor’s  case  [46]  is  shown  in  Fig.  16.  Using  the 
methodology  employed  in  the  previous  section  with  no  calibration  constants,  the  model  predicts  the  Sauter 
Mean  Diameter,  SMD,  with  reasonable  accuracy.  As  shown  in  Fig.  16,  there  is  a  steep  gradient  at  a  jet  speed 
around  U  «  20  m/s  for  the  Wu  et  al.  [52]  turbulent  jet  experiment.  Our  model  result  overlaps  with  that 
obtained  by  Hoyt  and  Taylor.  It  is  interesting  to  observe  that  Wu  et  al.’s  data  is  also  similar  to  our  result 
and  that  of  Hoyt  and  Taylor  at  about  U  «  20  m/s.  While  Wu  et  al.  noted  this  as  the  region  of  ‘uncertainty’, 
it  is  possible  that  the  rollup  motion  was  competing  with  the  turbulence  and  thus  the  perceptible  effect  of 
the  rollup  motion  appears  as  shown  in  Fig.  16. 

For  U  >  20  m/s,  differences  between  the  calculations  and  the  experiments  emerge.  It  is  known  that  linear 
analysis  [54]  overpredicts  the  droplet  size  (by  less  than  20%)  because  it  neglects  the  satellite  droplet  mass  due 
to  the  nonlinear  effect  [55,  56]  which  yields  the  multiple  crests  per  wavelength.  Ponstein’s  Eqn.  (41)  is  a  linear 
analysis  and  thus  it  also  tends  to  overpredict  the  droplet  size.  However,  a  20%  or  smaller  SMD  difference  does 
not  explain  the  difference  we  see  in  Fig.  16  at  a  higher  jet  speed.  This  is  due  to  the  fundamental  difference 
between  the  boundary  layer  instability  jet  and  the  turbulent  jet:  the  boundary  layer  instability  jet  is  scaled 
by  the  momentum  thickness  [42]  and  the  turbulent  jet  is  scaled  by  the  Kolmogrov  length  scale,  lk  or/and 
turbulence  eddy  characteristics  length  of  kinetic  energy,  h  [52],  Wu  et  al.  derived  the  empirical  formula 
using  the  ‘surface  kinetic  energy’  argument  which  gives  SMD  scaled  by  ~  1/U1A8.  Thus  the  governing 
length  scale  (i.e.  lk  and  h)  decrease  significantly  at  about  U  ~  20  m/s.  On  the  other  hand,  the  SMD  of 
the  boundary  layer  instability  jet  is  scaled  by  82 :  SMD  ~  l/U0'5  and  thus  its  change  with  respect  to  U  is 
relatively  moderate  as  shown  in  Fig.  16. 

Another  experiment  on  the  boundary  layer  instability  jet  is  available  by  McCarthy  and  Molloy  [57]. 
The  rollup  motion  causes  the  jet  to  be  atomized  when  the  circulation  is  large  enough  to  win  against  other 
competing  forces  such  as  viscous  or/and  capillary  forces.  In  Fig.  17-(a),  the  ‘stretching’  is  observed  due 
to  capillary  force  when  the  jet  is  atomized.  In  Fig.  17- (b),  the  computational  result  for  the  case  is  shown 
where  the  similar  ‘stretching’  of  the  capillary  force  is  observed.  For  the  high  speed  jet,  like  that  of  Hoyt 
and  Taylor  [46,  16],  the  effect  of  viscosity  and  capillary  force  is  of  little  importance  because  the  jet  is 
nearly  inviscid.  In  fact,  the  large  scale  motion  of  Hoyt  and  Taylor’s  jet  is  governed  by  the  Rayleigh  inviscid 
analysis  [54]  which  concludes  the  most  dominant  wavelength  to  be  A  =  4.51d. 

CONCLUSIONS 

A  fully  nonlinear  model  has  been  developed  to  simulate  primary  atomization  caused  by  boundary  layer 
instability  using  superposition  of  a  ring  vortex  with  a  potential  jet  flow.  The  axisymmetric  model  employs 
a  boundary  element  methodology  in  which  the  vorticity  in  the  boundary  layer  at  the  orifice  exit  is  used  to 
determine  ring  vortex  strength  and  radial  location  at  the  orifice  exit  plane.  Annular  ligaments  are  pinched 
off  the  surface  in  this  case;  a  secondary  linear  instability  analysis  due  to  Ponstein  is  used  to  predict  the 
fractionization  of  the  ligaments  into  individual  droplets. 

The  SMD  of  the  model  result  agrees  well  with  the  actual  droplet  size  of  Hoyt  and  Taylor’s  experiment. 
The  result  of  the  current  model  is  also  compared  with  the  experimental  data  of  Wu  et  al.  [52].  The  comparison 
confirms  that  the  current  model  predicts  the  droplet  size  satisfactorily. 
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Figure  1:  Typical  water  jet  into  air  in  the  atomization  regime.  Experimental  image  by  Hoyt  and  Taylor  [46]. 
Printed  under  the  permission  of  Journal  of  Fluid  Mechanics . 


Figure  2:  Schematic  of  the  initial  jet  geometry  indicating  computational  nodes  and  the  axisymmetric  ring 
vortex  at  the  orifice  exit  plane. 
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Figure  3:  Comparison  of  the  actual  physical  flow  condition  with  the  superposition  model  in  simulating 
boundary  layer  relaxation  downstream  of  the  orifice  exit  plane. 
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Figure  4:  A  9  x  9  matrix  decomposition  with  a  2  x  2  submatrix  for  a  2  x  3  process  grid 
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Figure  5:  Passing  the  solution  [A]  to  all  other  processors  is  carried  out  for  all  processors  to  share  the 
information  of  [X], 
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Figure  6:  Filter  transfer  function  vs.  wavenumber  for  various  filters 


Figure  7:  General  approximation  of  circulation  for  filament  vortex-ring  using  the  Kutta  condition  and 
vorticity-match  method. 


Figure  8:  (a)  Closeup  of  the  actual  Hoyt  and  Taylor’s  water  jet  [46].  Printed  under  the  permission  of  Journal 
of  Fluid  Mechanics,  (b)  The  closeup  of  the  model  result  for  Hoyt  and  Taylor’s  water  jet. 


Figure  9:  Application  of  Ponstein’s  [1]  theory  for  the  secondary  instability  of  a  pinch-off  vortex-ring 
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Figure  12:  Effect  of  the  initially  assumed  shape 
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ND  =  0. 1 971  +  0.7237  r.  +  0.8689 


-  Least  square  fit 

•  Model  result 


Figure  14:  Prediction  of  cicumferential  wave  number  (or  number  of  droplet)  due  to  circulation  around  the 
rotating  ring  pinched-off  from  the  main  liquid  stream. 


Figure  15:  Number  of  waves  vs.  core  thickness  of  pinch-off  vortex-ring 


U  [m/s] 


Figure  16:  Sauter  mean  diameter  comparison  for  various  jet  speed  U  in  [m/s]. 


(a)  Experiment 


(b)  Model  result 


Figure  17:  Liquid:  60%  glycerol  and  40%  water  by  weight,  p  =  103^,  p  =  llcP ,  a  =  0.0669^,  U  =  20  f , 
d  =  2.54 mm,  We  =  781,  Re  =  4750.  Comparison  between  experiment  [57]  and  model  results.  The  black 
circles  represent  the  location  and  the  relative  size  of  the  pinched-off  droplets.  Printed  under  the  permission 
of  Elsevier  Science. 
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6  Appendix  B  -  Jet  Stability  Theory 


Yoon*,  S.  S.  and  Heister,  S.  D.,  “Categorizing  Linear  Theories  for  Atomizing  Jets”. 
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Abstract 

This  paper  compares  and  contrasts  linear  stability  analyses  based  on  the  deformations  of  an  infinite  liquid 
column  and  due  to  boundary  layer  vorticity  imparted  to  the  free  surface  from  the  orifice  exit  plane.  The  bulk 
of  the  prior  works  which  date  back  to  Kelvin-Helmholtz  [1],  Rayleigh  [2],  Weber  [3],  Taylor  [4],  Ponstein  [5], 
Levich  [6],  Sterling  and  Sleicher  [7],  and  Reitz  and  Bracco  [8]  have  focused  on  the  liquid  column  analysis. 
Though  used  to  a  lesser  extent,  the  boundary  layer  instability  analysis  by  Brennen’s  [9]  can  also  be  used  to 
predict  the  dominant  wavelength  of  the  laminar  jet.  Differences  between  the  two  approaches  are  highlighted 
for  sample  injection  conditions  and  injector  geometries. 

1  Introduction 

The  linear  theory  associated  with  primary  atomization  of  a  liquid  jet  is  well  established  and  much  studied. 
More  recently,  nonlinear  analyses  have  served  to  compliment  the  linear  results  and  explain  the  presence  of 
satellite  droplets  which  have  frequently  been  observed  in  low-speed  jets.  At  higher  jet  speeds,  it  becomes 
more  difficult  to  assess  the  value  of  linear  theories  because  turbulence  and  secondary  atomization  become 
important  contributors  and  because  the  jet  surface  is  often  obscured  by  droplets.  In  spite  of  these  drawbacks, 
many  of  today’s  atomization  models  draw  heavily  from  the  linear  theories  based  on  a  column  of  liquid  exposed 
to  a  high-velocity  gas  shear  flow. 

The  biggest  drawback  from  these  theories  is  that  there  is  no  mechanism  to  include  effects  of  the  orifice 
geometry  on  the  possible  wavelengths  introduced  to  the  liquid  surface.  A  notable  exception  here  is  the  work 
of  Hoyt  and  Taylor  [10,  11,  12]  which  focused  on  the  boundary  layer  behavior  at  the  orifice  exit  plane  as  a 
significant  player  in  the  subsequent  instabilities  observed  in  the  jet.  Drawing  from  a  boundary  layer-based 
instability  analysis  due  to  Brennen  [9],  these  researchers  showed  that  the  wavelengths  observed  in  their 
carefully-controlled  experiment  were  in  agreement  with  those  predicted  using  Brennen’s  theory. 

The  aim  of  this  short  paper  is  to  compare  and  contrast  these  approaches  and  to  assess  potentially 
important  contributions  of  other  researchers  which  have  not  been  frequently  cited  in  prior  literature.  In  this 
context,  we  provide  a  brief  review  of  linear  stability  results  and  highlight  similarities  and  differences  in  the 
various  approaches.  Specific  examples  are  considered  to  illustrate  differences  between  liquid  column-based 
and  boundary  layer-based  approaches. 


2  Lord  Rayleigh’s  Analysis 


Probably  the  most  famous  and  poineering  stability  analysis  for  the  liquid  jet  was  developed  by  Lord 
Rayleigh  [2].  Rayleigh  considered  the  infinitely  long  inviscid  column  of  liquid  with  negligible  influence 
from  the  gas  phase.  He  hypothesized  the  infinitesimal  disturbance  will  cause  the  jet  to  breakup  under  a 
capillary-based  instability.  The  famous  dispersion  relation  he  obtained  is  as  follows: 
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where  w  =  wr  +  iwt  (i.e.  wr=growth  rate,  i  =  V^T,  and  u>i=frequency  of  oscillation),  a=surface  tension 
of  the  liquid,  p;=liquid  density,  a=orifice  radius,  k= wave  number=27r/A  (i.e.  A= wavelength),  I\ {ka)  and 
I0(ka)  are  modified  Bessel  functions  of  the  first  kind.  By  expanding  the  Bessel  functions  in  a  power  series 
and  computing  the  maximum  of  the  w  vs.  ka  curve,  he  obtained  the  result: 
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for  which  the  corresponding  wavnumber  and  wavelength  are: 

ka  =  0.696  «  0.7  \  =  4.51  d 

There  are  numerous  experimental  confirmations  of  Rayleigh’s  wavelength  appearing  in  low-speed  jets.  Many 
researchers  have  extended  the  result  to  the  fully  nonlinear  regime  and  assumed  that  the  jet  will  actually 
fragment  into  sections  4.51d  in  length  which  leads  to  a  predicted  droplet  diameter  of  about  1.89d  as  shown 
in  Fig.  1.  While  such  conditions  can  be  generated  with  a  carefully  controlled  perturbation  in  an  experiment, 
more  commonly  the  jet  is  shown  to  bifurcate  into  two  drops  per  wavelength;  the  presence  of  the  satellite 
drops  are  not  predicted  from  linear  theory,  but  numerous  nonlinear  results  [13,  14, 15,  16,  17]  have  confirmed 
their  presence  and  agree  well  with  size  measurements  from  a  variety  of  experiments. 

The  Rayleigh  jet  is  of  course  the  simplest  of  all  cases  in  that  aerodynamic  interactions  with  the  gas  are 
neglected.  The  low  jet  velocities  associated  with  this  flow  regime  also  imply  that  instabilities  emanating 
from  the  boundary  layer  inside  the  orifice  are  necessarily  small.  For  these  reasons,  there  is  good  agreement 
between  experiment  and  theory  for  a  variety  of  orifice  designs. 


3  Weber’s  Equation 


Weber  [3]  extended  Rayleigh’s  analysis  by  adding  the  effect  of  viscosity  of  the  jet  which  gives: 

^  +  ^{ka)2w  =  2^(1  “  W 

The  coefficient  of  the  w  term  accounts  for  viscous  effects.  McCarthy  and  Molloy  [18]  provided  the  maximum 
growth  rate  of  Weber’s  Eq.  (4). 
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The  corresponding  most  dominant  wavelength  is 
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Nonlinear  simulations  for  viscous  [19]  and  inviscid  [13,  20]  flows  show  comparable  results  for  low-speed  jets 
as  well;  the  presence  of  viscosity  only  seems  to  decrease  the  rate  at  which  the  instability  grows,  not  the  shape 
of  the  wave. 


4  Sterling-Sleicher  Equation 


The  greatest  shortcoming  of  Rayleigh  and  Weber  equations  is  due  to  the  neglect  of  aerodynamic  forces  on 
the  liquid  jet.  Sterling-Sleicher  [7]  (hereafter  referred  as  ‘SS’)  had  developed  the  dispersion  equation  that 
takes  into  accout  the  aerodynamic  effects;  the  general  SS  equation  is: 
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In  the  absence  of  gas  phase  and  viscosity  (i.e.  e  =  p  =  0),  the  SS  Eq.  (7)  reduces  to  Rayleigh’s  result.  In 
addition,  the  general  SS  Eq.  (7)  can  be  written  as  the  following  for  the  inviscid  case: 
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Sterling-Sleicher  considered  the  case  where  £  =  ka  <  1.0  which  implies  the  wavenumber  is  small  or  the 
wavelength  (i.e.  A  =  2tt /k)  is  large  enough  that  it  is  on  the  order  of  the  nozzle  radius,  a.  In  this  case,  the 
general  SS  equation  is  simplified  substantially: 
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It  should  be  noted  that  a  typoghraphical  error  is  present  in  the  original  Sterling-Sleicher  paper  in  the  ‘e’ 
term.  It  should  be  U 2  while  it  was  typed  as  just  U.  Here  U  is  the  ‘constant’  and  ‘uniform’  jet  speed.  In  the 
limit,  of  £  — >  oo  for  an  inviscid  case,  the  SS  Eq.  (7)  is  essentially  the  same  as  the  Kelvin-Helmholtz  equation, 
which  is  used  to  predict  the  smaller  wavelengths  (detailed  proof  of  this  claim  is  available  in  the  following 
section). 


5  Reitz-Bracco  Equation 

Further  extention  of  the  SS  equation  was  presented  by  Reitz-Bracco  [21]  with  the  assumption  that  the 
liquid  jet  velocity  is  also  a  function  of  the  radial  direction  (i.e.  U  =  Ug(r)),  as  in  the  Orr-Sommerfeld 
equation  [22,  23].  This  equation  of  Reitz-Bracco  (hereafter  referred  as  ‘RB’)  seems  to  be  the  most  general 
dispersion  expression  available  in  the  literature  for  the  axisymmetric  case. 
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It  should  be  noted  that  Levich  [6]  had  performed  a  similar  analysis  the  result  of  which  is  exactly  the  same 
as  the  RB  equation  except  that  Levich  had  omitted  the  terms  pgw2  and  pgikUw  of  gas  by  assuming  the 
two  terms  are  negligible  compared  to  the  terms  piw2  and  piikUw  of  liquid.  In  the  absence  of  viscosity  and 
the  gas  phase  (i.e.  v  =  e  =  0),  the  RB  Eq.  (11)  also  recovers  Rayleigh’s  result.  Further,  the  RB  Eq.  (11) 
becomes  the  following  if  inviscid  (i.e.  v  =  0): 
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This  is  exactly  the  same  as  the  inviscid  case  of  the  SS  Eq.  (9)  except  for  the  coefficient  of  in2.  This  difference 
originated  from  Sterling-Sleicher’s  uniform  jet  speed  (i.e.  U  =  constant )  assumption  which  differs  from 
Reitz-Bracco’s  jet  velocity  approximation  with  the  gas  velocity  as  a  function  of  the  radial  direction.  It  is 
interesting  to  note  that  both  the  inviscid  SS  Eq.  (9)  and  that  of  RB  Eq.  (13)  are  reduced  to  the  following 
expression  in  the  limit  as  £  — ►  oo: 

w2  =  eU2k2  -  —  (14) 

Pi 

where  we  have  used  the  approximations  lim^_oo  Io(0/I i(£)  =  0,  lim^_oo  Ko(Q/K i(£)  =  1  (see  Pearson  [24]), 
and  pi  pg.  This  Eq.  (14)  is  essentially  the  Kelvin-Helmholtz  equation  [25,  1]. 
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6  Ponstein’s  Equation 

In  this  section,  the  linear  theory  of  Ponstein  [5]  is  discussed.  His  work  published  in  1959  in  Applied  Scientific 
Research,  Deutch  Journal,  has  not  been  well  recognized  in  the  atomization  community.  The  work  not  only 
had  extended  Rayleigh  s  [2,  26]  analysis  to  include  gas-phase  effects,  but  also  considered  column  rotation 
(swirl)  in  an  analysis  published  long  before  the  Sterling-Sleicher  (1975),  and  even  before  Levich  (1962). 

Ponstein  had  considered  two  cases:  a  rotating  liquid  column  in  gas-phase  and  a  rotating  bubble  (or  gas) 
column  in  a  liquid  surrounding  for  the  second  case.  A  uniform  liquid  column  in  vacuum  is  well  known 
by  Rayleigh  [2]  who  predicted  the  most  dominant  wavelength,  A  =  4.5  Id.  Rayleigh  [26]  also  considered  a 
uniform  bubble  column  in  liquid  whose  solution  is: 


u2  = 


pia- 


(1  —  k2a2)£ 


Ko(Z) 


(15) 


This  equation  predicts  a  most  unstable  wavelength,  A  =  6.48d.  For  an  axisymmetric  rotating  bubble  column 
(based  on  eut),  Ponstein  gives  the  following  result: 


LJ 
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(16) 


where  T  is  the  circulation  around  the  ring  (or  column)  which  can  be  estimated  as  T  =  (2tt a)VB  from 
Saffman  [27],  Here  Vg  is  the  tangential  velocity  of  the  ring  surface.  For  a  non-rotating  case  (i.e.  T  =  0), 
Eq.  (16)  recovers  the  Rayleigh’s  result  in  Eq.  (15).  In  this  case,  circulation  has  a  stabilizing  influence  as 
indicated  by  the  negative  sign  on  the  T  term.  The  faster  it  rotates,  the  more  stable  the  bubble  ring  is. 
The  detailed  discussion  of  Eq.  (16)  is  available  in  Lundgren  and  Mansour  [28]  where  they  had  modeled  the 
evolution  of  the  bubble  vortex-ring  using  boundary  integral  method.  An  interesting  example  of  the  bubble 
vortex-ring  of  Dolphin  is  discussed  in  Shariff  [29]. 

Ponstein  gives  the  following  result  for  the  second  case  he  had  considered,  a  rotating  liquid  column  in  gas: 


u2  = 


— 3  (1  -  k2a2)  +  (1  -  e) 


+  eU2k2 


(17) 


If  we  consider  non-rotating  (i.e.  T  =  0)  and  non-aerodynamic  effect  (i.e.  U  =  e  =  0),  Rayleigh’s  result  is 
recovered.  Here,  circulation  has  a  destabilizing  effect  as  indicated  by  the  positive  sign  on  the  T  term.  The 
faster  the  column  rotates,  the  more  unstable  it  becomes.  Increasing  gas  density  e  serves  to  aid  in  stabilizing 
the  column  circulation  term,  but  destabilizes  the  dominant  aerodynamic  ( U 2)  term. 

Considering  the  non-rotating  case  with  aerodynamic  effect,  Ponstein’s  Eq.  (17)  can  be  written  as: 


For  £  <  1.0,  it  is  known  that  Ji(f )//<,(£)  «  (£)/2.  Applying  this  identity,  Eq.  (18)  is  re-written  as: 
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This  result  is  exactly  the  same  as  the  inviscid  case  of  the  dispersion  relation  derived  by  Sterling-Sleicher  in 
Eq.  (10). 


7  Taylor’s  Equation 

Taylor  [4]  derived  the  following  relation  based  on  his  aerodynamic  theory: 

m  =  ]f^2x9{x)  <2°) 


43 


where  x  is  the  nondimensional  wavelength  of  Reitz  and  Bracco  [21]  and  g  —  g{x)  is  the  function  of  the 
Taylor’s  number,  Ta. 


ak 


(21) 


with  Re  =  piUd/g,  and  Weltd  =  piU2d/a.  The  RB  Eq.  (11)  in  the  limit  of  £  ->  oo  with  viscosity  (i.e.  v^O) 

is  as  follows:  3  - 

(w  +  2 vk2)2  +  —  -  4 u2k\k2  +  -  +  +  iUkf  =  0  (22) 

V  '  Pi  V  v  pi 

This  equation  is  modified  in  the  form  of  the  Taylor’s  Eq.  (20)  assuming  pj  >  pg, 


w 

Uk 


(23) 


where  _ 

AM  =  x-2(2^v^7f-2rf»-2^v)  (24) 

If  we  hypothesized  the  RB  Eq.  (23)  were  the  same  as  Taylor’s  Eq.  (20),  the  following  relation  can  be  set: 

2  x  g(x)  =  h(x)  (25) 


Reitz  and  Bracco  [21]  also  defined  a  function  f(x)  =  x  g(x).  The  numerical  solution  of  Eq.  (23)  was  given 
by  Reitz  and  Bracco.  The  most  dominant  nondimensional  wavelength  xm  was  found  to  be  1.5  and  thus 
f(xm)  -  v^S/6  «  0.2887  as  Ta  >  1.0  (see  Bracco  [30]).  Though  the  value  of  f(xm)  =  \/3/6  was  first  found 
by  G.  I.  Taylor  [4]  and  was  mentioned  in  Ranz’s  paper  [31],  it  seems  that  the  most  dominant  nondimensional 
wavenumber  xm  =  1.5,  which  gives  f(x  m)  —  \/3/6,  as  noted  by  Reitz  and  Bracco  [21].  The  same  results 
have  been  found  using  KH  Eq.  (14)  which  is  significantly  simpler  than  Eq.  (23).  KH  Eq.  (14)  can  be  written 
in  the  form  of  Taylor’s  Eq.  (20): 


Note  that  KH  Eq.  (14)  is  valid  under  inviscid  (i.e.  v  =  0)  assumption.  If  this  were  the  same  as  Taylor’s 
Eq.  (20),  2 xg{x)  =  (1  -  £)1/2  can  be  set: 


(27) 


The  f(x)  and  g(x)  as  a  function  of  the  nondimensional  wavelength  x  are  plotted  in  Fig.  2.  The  g{xm)  can 
be  found  by  taking  derivative  of  Eq.  (27)  with  respect  to  x. 


dg(x)  _  1  —  2xy 
dx  Axzy1/2 


(28) 


where  y=l-~.  Let  1  -  2xy  =  0  to  find  the  most  dominant  nondimensional  wavenumber,  xm.  This  gives 
the  following  analytical  results: 

\/3 

Xm  =  1.5  f(xm)  =  Xm  g{xm)  =  -g"  (29) 

This  shows  that  the  viscosity  is  of  little  importance  since  no  difference  is  seen  in  xm  and  f{xm)  values 
between  Reitz  and  Bracco’s  numerical  result  and  our  analytical  results  using  KH  Eq.  (14)  (see  Fig.  2). 

Wu  et  al  [8]  assumed  that  the  initial  droplet  diameter,  D,  might  be  proportional  to  the  length  of  the 
most  unstable  wavelength,  A 

D  —  BX  (30) 
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B  is  a  constant  that  is  empirically  approximated.  Wu,  Reitz  and  Bracco  [8]  gave  B  «  4.5,  based  on  the 
curve-fit  of  their  experimental  data.  Applying  this  Eq.  (31)  to  Hoyt  and  Taylor’s  case  [11]  which  is  discussed 
in  detail  in  the  next  section  (i.e.  aHio  =  0.0734 kg/s2,  pair  =  1.23 kg/m3,  U  =  21m/s,  and  d  =  6.35mm, 
thus  Weg  =  46.9),  the  most  unstable  wavelength  is  A  w  d/5.  The  wavelength  overpredicts  the  value  observed 
in  Hoyt  and  Taylor’s  experiment  (i.e.  A  «  d/13.8).  It  should  be  noted  that  Eq.  (31)  is  capable  of  producing 
very  small  wavelengths  with  large  U  and  pg,  which  is  the  case  of  Wu  et  al  [8];  this  aerodynamic  theory  works 
well  for  their  atomizing  jets  of  the  shear  layer  driven  instability. 


8  Boundary  Layer  Instability  Analysis 


None  of  the  conventional  linear  stability  theories  can  account  for  the  largely  axisymmetric  waves  observed  by 
Hoyt  and  Taylor  [11]  in  their  famous  1977  experiments.  In  this  experiment,  a  carefully  machined  nozzle  was 
used  to  provide  a  favorable  pressure  gradient  to  a  reasonably  high-speed  jet  such  that  turbulence  effects  were 
minimized.  The  resultant  waves  formed  on  the  surface  of  the  water  jet  showed  wavelengths  A  «  d/13.8  which 
is  not  predicted  well  by  the  linear  theories  discussed  above.  These  researchers  theorized  that  the  thickness  of 
the  boundary  layer  at  the  orifice  exit  plane  could  play  a  significant  role  in  explaining  the  observed  wavelength. 
It  was  suspected  that  the  point  of  inflection  due  to  change  in  velocity  profile  from  a  no-slip  to  a  free-surface 
edge  condition  was  responsible  for  the  instability  and  the  subsequent  wave  growth.  Shkadov  [32]  was  another 
researcher  working  on  this  theory  at  the  time. 

For  this  reason,  we  review  the  classical  boundary  layer  instability  based  on  the  Orr-Sommerfeld  equation. 
Using  (i)  the  solution  by  Betchov  and  Criminale  [33]  for  the  fully  developed  2D  steady  parallel  laminar  flow, 
(**)  the  perturbation  method  [22],  (Hi)  and  the  linearization,  the  2D  Navier-Stokes  equations  can  be  written 
as  follows: 


v  -  \v" 

1  .  ,, 

H — ri >u  = 
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25"  -  a2v  -  -v"" 

a2 

a 2 

Res2  io. 

a 

where  u  =  u(y)  is  the  axial  velocity  profile  in  y  direction,  c  =  cr-Hct  is  the  growth-rate,  a  is  the  wavenumber, 
v  =  v(y)  is  the  amplitude  of  the  radial  perturbation  in  y  direction,  and  Re&2  =  Ud2/v.  This  is  called  ‘Orr- 
Sommerfeld  equation.  Sometimes  different  notations  (i.e.  v(y)  =  —ict<f>( y))  are  used  by  researchers  (see 
Panton  [23]  and/or  Schlichting  [22]).  It  should  also  be  noted  that  <j>(y )  is  not  velocity  potential  but  complex 
amplitude  function  for  the  streamfunction.  There  is  no  analytic  solution  to  Eq.  (32)  but  the  numerical 
solution  for  a  wide  range  of  values  of  frequency  and  ReS2  is  available  by  Jordinson  [34], 

Considering  the  inviscid  case  of  Orr-Sommerfeld  Eq.  (32)  (i.e.  Re$2  — >  oo),  the  following  equation  is 
obtained: 


(u-c) 
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This  is  called  ‘Rayleigh  Equation’.  Lord  Rayleigh  (1880-1913)  provided  a  theorem  based  on  his  Eq.  (33), 
“Velocity  profiles  with  points  of  inflection  are  unstable” .  Rayleigh  proved  that  the  presence  of  a  point-of- 
infection  is  a  necessary  (though  it  is  not  sufficient)  condition  for  the  appearance  of  unstable  waves.  This  is  a 
powerful  statement  as  it  can  be  used  to  classify  the  flow  regime:  laminar  to  turbulent  flows.  For  example,  the 
velocity  profile  contains  no  point-of-inflection  under  favorable  pressure  gradient  (i.e.  dP/dx  <  0).  Generally, 
it  is  fair  to  state  that  the  flow  is  laminar  in  such  case.  On  the  other  hand,  point-of-inflection  can  be  found 
in  adverse  pressure  gradient  (i.e.  dP / dx  >  0)  where  the  flow  is  sometimes  unstable  and  eventually  this  may 
lead  to  turbulent  flow. 

The  influence  of  gas  density,  and  hence  aerodynamic  forces  on  the  interface,  has  been  widely  investigated 
experimentally.  Reitz  and  Bracco  [21]  observed  a  substantial  difference  in  the  atomization  mechanism  when 
the  liquid  jet  was  injected  in  different  gases  (i.e.  pNi  =  6  kg/m3  and  pXe  =  23  kg/m3).  Wu  et  al  [35] 
have  reported  a  change  in  droplet  size  for  primary  atomization  when  a  different  gas  density  was  tested  for 
the  same  liquid  jet.  From  this  basis  one  may  conclude  that  the  aerodynamic  interaction  at  the  surface  of 
the  jet  does  alter  wave  growth,  i.e.  that  aerodynamic  forces  are  of  sufficient  magnitude  to  contribute  to 
the  instability.  Shkadov’s  theory  [32]  states  that  the  velocity  profile  of  liquid  is  independent  of  gas  density, 
which  seems  contradictory  to  the  experimentally  observed  trend.  This  dilemma  is  probably  due  to  the  high 
gas-liquid  density  ratio  which  causes  the  growth  rate  due  to  the  aerodynamic  term  to  dominate  as  compared 
to  the  growth  rate  due  to  the  vortices,  thereby  invalidating  Shkadov’s  theory.  In  addition,  the  droplet  sizes 
are  strongly  dependent  on  secondary  atomization  at  high  jet  speed  and  high  gas  density  conditions. 
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Clearly  the  influence  of  the  gas  is  regime  dependent,  but  the  fundamental  instability  mechanism  is  present 
in  all  regimes.  The  fundamental  questions  which  should  be  asked  are,  “Will  instability  occur  in  the  absence 
of  gas  phase?  If  so,  where  is  that  instability  mechanism  originating  from?” 

Hoyt  and  Taylor[10]  did  not  observe  differences  in  wave  structure  with  differing  external  airflows  in  their 
experiments.  The  axisymmetrically  disturbed  short  wavelength  observed  near  the  nozzle  exit  is  present 
regardless  of  magnitude  or  direction  of  the  air  velocity.  Hoyt  and  Taylor  concluded  that  this  phenomenon 
has  “no  discernible  effect  of  relative  air  velocity”.  In  spite  of  Hoyt  and  Taylor’s  effort  in  this  regard,  their 
result  has  not  well  been  recognized  in  the  atomization  community.  For  many  years,  the  notion  that  vorticity 
at  the  nozzle  exit  being  responsible  for  the  atomization  has  been  over-shadowed  by  the  aerodynamic  linear 
theory  and  the  experiment  which  showed  a  significant  difference  in  spray  structure  with  the  change  in  gas 
density. 

Shkadov  [32]  provides  the  solution  of  the  Rayleigh’s  Eq.  (33)  and  proves  that  the  amplitude  of  surface 
waves  grows  in  the  downstream  direction,  as  the  jet  velocity  profiles  relaxes.  This  eigenvalue  problem  is  also 
solved  by  Brennen  [9].  Brennen  provided  the  boundary  layer  instability  analysis.  He  considered  separated 
boundary  layer  flow  over  the  planar  plate  using  Gaussian  velocity  profile.  This  resulted  in: 

7  =  27r/ ”  (34) 


where  7  is  the  nondimensional  frequency,  /  is  the  dimensional  frequency  in  [Hz],  S2  is  the  momentum 
thickness  in  [m],  and  U  is  the  speed  of  the  uniform  flow  in  [m/s].  Brennen  concluded  that  7  =  0.175  was  to 
be  the  nondimensional  frequency  which  would  give  maximum  amplification  at  the  flow  separation  point. 
Applying  7  =  0.175  to  Hoyt  and  Taylor  case  [11] 


7  TJ  _  ( 0.175  \  /  21  m/s  \ 

27t52  \  27r  1.2x10 ~5m  J 


48, 741 Hz 


(35) 


Note  that  the  boundary  layer  momentum  thickness,  <52,  can  be  approximated  using  Blasius  solution  [36]  for 
laminar  flow  assuming  that  S2  «  a: 


Si  _  1.721  62  _  0-664 

X  yj ~Rg%  r  y/ 


(36) 


where  Rex  —  and  <5i  is  the  displacement  thickness.  Assuming  that  the  wave  speed  is  about  the  same  as 
that  of  the  liquid  jet  U, 

U  21  m/s  d  ,  . 

Aw  —  = - - —  =  0.43  mm  =  — —  (37 j 

/  48741ifz  14.8 


While  this  was  a  theoretically  predicted  value,  Hoyt  and  Taylor’s  experimental  observation  of  the  axisym¬ 
metrically  disturbed  wavelength  was  about  A  w  0.46  mm  =  d/13.8  as  shown  in  Fig.  3  (note:  d  =  6.35  mm). 
The  comparison  between  theory  and  experiment  was  excellent. 

McCarthy  and  Molloy  [18]  ’s  experiment  provides  further  confirmation  of  the  role  of  the  momentum 
thickness  in  wave  formation  at  the  orifice  exit  plane.  They  varied  the  nozzle-to-diameter  ratio  (i.e.  I/d) 
and  observed  the  effect  of  l/d  on  ‘laminar’  jet  structure.  In  their  paper  [18],  it  was  not  mentioned  that  the 
jets  had  distinctive  axisymmetric  waves  for  different  orifice  designs,  yet  the  experimental  images  do  contain 
these  structures  as  shown  in  Fig.  4  excerpted  from  their  paper.  In  the  further  downstream,  the  jets  become 
turbulent  (see  Case  l/d  =  5  and  10  in  Fig.  5)  and  result  in  the  primary  atomization. 

To  compare  the  McCarthy  and  Molloy  results  with  Brennan’s  theory,  we  assume  boundary  layer  devel¬ 
opment  inside  the  passage  can  be  approximated  by  boundary  layer  growth  on  a  flat  plate.  Clearly,  this 
assumption  becomes  poor  when  the  boundary  layer  is  a  substantial  fraction  of  the  orifice  radius,  but  one 
could  use  a  numerical  analysis  or  more  elaborate  theory  to  more  accurately  ascertain  momentum  thicknesses 
at  the  orifice  exit  plane.  The  errors  result  from  freestream  pressure  gradients  which  develop  as  the  boundary 
layer  builds  in  the  passage.  A  thorough  review  of  turbulent  momentum  thickness  is  available  by  Klein  [37]. 
Given  this  caveat,  the  predicted  6 2  values  using  this  technique,  when  implemented  in  Brennen’s  equation, 
show  good  agreement  with  the  observed  wavelengths  from  the  McCarthy  and  Molloy  experiments  as  pre¬ 
sented  in  Table  1.  One  can  also  note  from  Table  1  that  the  S2  values  are  all  quite  small  relative  to  the  orifice 
radius,  thereby  lending  credence  to  the  simple  approach  of  assuming  flat  plate  boundary  layer  growth. 
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Table  1:  Summary  for  McCarthy  and  Molloy’s  Experiment  [18] 


l/d * 

Si/d 

/"  \Hz] 

A« 

1 

1/40.04 

1/103.8 

46 

22758 

d/2.9 

d/2.8 

5 

23738 

1/17.90 

1/46.41 

102 

10178 

d/1.3 

d/1.4 

10 

47477 

1/12.66 

1/32.82 

145 

7197 

1.1  d 

1.0  d 

*  l/d= nozzle  length  to  diameter  ratio 


t  Rex  =  VI jv  M  Brennen’s  Eqs.  (35)  and  (37),  respectively 
note:  Blasius  solution  is  used  for  6 1  and  62  estimation 


9  Comparison  of  the  Linear  Models 

It  is  useful  to  compare  the  wavelengths  predicted  by  the  liquid  column  theories  and  the  boundary  layer 
theory  for  the  same  atomizer.  In  Fig.  6,  wavelength  predicted  by  Brennen’s  theory  is  plotted  as  a  function 
of  jet  speed  for  various  l/d  and  is  compared  with  the  SS  and  KH  liquid  column-based  analyses. 

The  results  obtained  by  Sterling-Sleicher  Eq.  (10)  and  that  obtained  by  Kelvin-Helmholtz  Eq.  (14)  are 
essentially  the  same  at  the  higher  jet  speeds  indicating  that  the  KH  mechanism  is  dominant  for  atomizing 
jets.  The  Reitz  and  Bracco  Eq.  (13)  (though  not  plotted  in  Fig.  6)  gives  results  very  similar  to  the  Sterling 
and  Sleicher  result  as  noted  in  prior  discussion.  These  column-based  results  show  a  strong  dependence  on  jet 
velocity  at  the  lower  speeds,  with  an  asymptotic  behavior  at  the  high  jet  speeds.  The  boundary  layer-based 
results  are  also  provided  in  Fig.  6  for  various  orifice  lengths.  These  results  show  a  more  modest  variation 
in  critical  wavelength  at  the  lower  injection  velocities  with  a  similar  asymptotic  behavior  at  the  high  jet 
speeds.  Results  from  the  boundary  layer-based  analysis  show  a  strong  influence  of  orifice  Ijd,  while  the 
column-based  analysis  does  not  include  this  parameter.  In  general,  the  boundary  layer-based  results  predict 
smaller  wavelengths  at  the  low  jet  speeds  and  larger  wavelengths  at  the  higher  jet  speeds  as  compared  to 
the  column-based  results. 

Unfortunately,  the  influence  of  turbulence,  especially  at  the  higher  jet  speeds,  makes  it  difficult  to  perform 
experiments  to  compare  the  two  theories.  The  axisymmetrically  perturbed  waves  are  not  observable  when 
the  jet  flow  is  already  turbulent  flow  at  the  nozzle  exit  (see  Figure- 1  of  Wu  et  al  [38]).  However,  the  flow 
at  Hoyt  and  Taylor  nozzle  exit  is  nearly  laminar  flow  as  the  flow  is  laminarized  through  highly  contracted 
nozzle  geometry  under  the  favorable  pressure  gradient  [11].  Their  carefully  conducted  experiment,  which 
had  substantially  reduced  all  other  possible  perturbation  sources  (i.e.  such  as  turbulence  and  cavitation), 
has  made  the  observation  of  the  boundary  layer  waves  possible.  For  their  water  jet  (laminar  at  the  nozzle 
exit),  the  Reynolds  number  based  on  the  momentum  thickness  (ReS2)  is  about  225  while  the  Blasius  solution 
indicates  that  the  flow  should  be  turbulent  for  Re$2  >  200. 

This  dilemma  is  probably  because  that  the  favorable  pressure  gradient  has  shifted  the  shape  factor  (i.e. 
H  =  Si/S2)  and  thus  the  critical  Reynolds  number  (Jte*2fCrit)  has  increased.  It  should  be  noted  that  a  slight 
decrease  in  H  can  result  in  a  substantial  increase  in  Res2>Crit  [36,  39].  At  this  stage  of  the  research,  it  is 
not  clear  where  the  Res2}Crit  lies  for  the  liquid  jet.  In  Fig.  7,  we  have  considered  a  laminar  water  jet  into 
an  air.  It  is  shown  that  Re&2  increases  with  increasing  l/d.  In  reality,  it  is  not  guaranteed  that  the  flow  at 
the  nozzle  exit  is  laminar  for  high  l/d.  However,  one  will  have  a  better  chance  of  observing  the  boundary 
layer  waves  if  Re&2  <  Res2iCrit •  It  certainly  would  be  interesting  to  generate  additional  experimental  data 
for  these  conditions. 


10  Conclusions 

Linear  analyses  based  on  liquid  column  and  boundary  layer-based  methodologies  are  compared  and  con¬ 
trasted.  The  boundary  layer-based  analysis  has  the  advantage  that  it  addresses  the  orifice  geometry  as  a 
primary  influence  in  determining  the  most  unstable  wavelength.  The  column-based  theories  are  shown  to 
be  largely  equivalent  in  the  atomization  regime  due  to  the  dominance  of  the  Kelvin-Helmholtz  instability 
mechanism  at  these  conditions.  In  performing  a  sample  comparison  of  the  two  methods  for  a  water  jet,  the 
boundary  layer-based  scheme  tends  to  predict  smaller  critical  wavelengths  at  low  speed  conditions  and  larger 
critical  wavelengths  at  high  injection  speeds.  While  the  column-based  analysis  has  received  more  attention 
in  the  community,  the  boundary  layer-based  approach  is  shown  to  have  merit  based  on  limited  experiments 
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aimed  at  addressing  this  mechanism. 
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Figure  1:  Typical  low  speed  jet  in  the  Rayleigh  regime.  Experimental  image  by  Adam  Hart 
Printed  under  the  permission  of  Adam  Hart-Davis. 


Figure  2:  The  functions  g(x)  by  Taylor  [4]  and  f{x)  =  xg(x)  by  Reitz  and  Bracco  [21]  as  a  function  of 
nondimensional  wavenumber  x  =  (paU2)/(ak). 


Figure  3:  Typical  water  jet  into  air  in  the  atomization  regime.  Experimental  image  by  Hoyt  and  Taylor  [10]. 
which  shows  the  most  dominant  wavelength  A  =  d/13.8  while  Brennen’s  theory  predicts  \B  =  d/14.8. 
Printed  under  the  permission  of  Journal  of  Fluid  Mechanics. 
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Figure  4:  McCarthy  and  Molloy’s  experiment  [18]  for  l/d=  10,  5,  and  1:  The  most  dominant  wavelength 
appears  subsequent  to  laminar  region  which  can  be  scaled  by  the  Brennen’s  [9]  theory.  Printed  under  the 
permission  of  Elsevier  Science. 
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Figure  5:  McCarthy  and  Molloy’s  experiment  [18]  for  l/d=  10,  5,  and  1  in  the  further  downstream.  Droplet 
size  is  not  scaled  by  the  most  dominant  wavelength.  Printed  under  the  permission  of  Elsevier  Science. 


U  [m/s] 

Figure  6:  Variation  of  the  theoretically  predicted  boundary  layer  waves  as  well  as  that  of  Sterling-Sleicher 
Eqs.  (10)  and  Kelvin-Helmholtz  Eq.  (14)  as  a  function  of  jet  speed  for  various  l/d.  Water  jet  into  air  is 
considered:  v  =  1.12  x  10~6m2/s  and  d  =  l  mm.  2nd  wind-induced  regime  for  26.76  m/s  <U  <  48.86  m/s 
and  atomization  regime  for  U  >  48.86  m/ s. 


7  Appendix  C  -  Acoustic  Interactions  with  Liquid  Jet 


Yoon*,  S.  S.  and  Heister,  S.  D.,  “Modeling  Atomizing  Jet  due  to  Boundary  Layer  Insta¬ 
bilities”  . 
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Abstract 

A  nonlinear  model  has  been  developed  to  assess  the  time  evolution  of  an  axisymmetric  liquid  jet  using  a 
boundary-element  method.  Vorticity  transported  from  the  boundary  layer  in  the  orifice  passage  to  the  free 
surface  is  modeled  using  a  potential  ring  vortex  placed  at  the  orifice  exit  plane.  The  vortex  strength  is 
uniquely  determined  from  the  Kutta  condition  and  information  regarding  the  boundary  layer  thickness  at 
the  orifice  exit  plane.  It  is  shown  that  primary  breakup  can  occur  even  in  the  absence  of  the  gas  phase. 
Using  a  secondary  stability  analysis  after  Ponstein  [1],  the  size  of  the  droplets  is  estimated  based  on  the  size 
of  the  ring-type  structures  shed  from  the  periphery  of  the  jet.  Computed  droplet  sizes  are  in  reasonable 
agreement  with  experimental  data,  although  turbulence  effects  obscure  some  comparisons. 

1  Introduction 

The  hydrodynamic  instability  of  liquid  jet  flow  has  been  of  great  interest  to  fluid  dynamicists  for  more  than 
a  hundred  years.  Despite  the  numerous  analytical,  experimental  and  numerical  studies  over  the  years,  the 
jet  flow  is  too  complex  to  be  understood  completely.  While  the  flow  in  the  first  few  diameters  downstream 
from  the  nozzle  exit  greatly  depends  on  the  internal  flow  effect,  the  relative  motion  between  the  main  liquid 
stream  and  the  gas  plays  an  important  role  farther  downstream.  Droplets  pinched  off  from  the  liquid  stream 
experience  drag  and  therefore  result  in  breakup.  In  addition,  collision  between  one  droplet  and  another 
enhances  the  complexities  of  the  jet  atomization  process.  These  phenomena  are  greatly  dependent  on  the 
initial  flow  condition  near  the  nozzle  exit. 

Two  upstream  flow  conditions  that  may  affect  the  initial  flow  condition  at  the  nozzle  exit  are  turbulence 
and  cavitation.  While  DeJuhasz  [2]  claimed  that  turbulence  may  be  the  most  important  factor  in  jet  breakup 
process,  it  was  later  shown  by  Bergwerk  [3]  that  the  turbulence  eddy  viscosity  in  the  applicable  range  of 
Reynolds  numbers  is  not  large  enough  to  cause  the  disintegration  of  the  jet.  Bergwerk  suspected  that 
cavitation  was  the  main  source  that  produces  an  amplitude  large  enough  to  cause  the  jet  breakup. 

In  the  absence  of  cavitation  and  the  substantial  reduction  of  turbulence  fluctuations  through  the  use 
of  a  nozzle  geometry  promoting  highly  favorable  pressure  gradients  [4]  atomization  is  still  known  to  occur. 
This  suggests  that  there  are  other  mechanisms  that  lead  to  the  disintegration  of  the  jet.  Rupe  [5]  observed 
the  velocity  profile  relaxation  has  a  key  role  in  influencing  the  jet  breakup.  The  boundary  layer  instability 
analysis  by  Shkadov  [6]  predicted  the  unstable  short  wavelength  of  the  free  surface  of  the  jet.  In  the  sudden 
absence  of  nozzle  wall  and  ignoring  jet  divergence  (i.e.  parallel  flow),  the  jet  velocity  at  the  free  surface 
will  start  to  accelerate  from  the  stagnant  condition  consistent  with  the  no  slip  boundary  condition  inside 
the  nozzle.  This  condition  results  in  a  point  of  inflection  in  the  velocity  profile,  which  is  inviscidly  unstable 
according  to  Rayleigh’s  theorem  (see  Schlichting  [7]  and  Panton  [8]);  vortices  start  to  form  at  the  point  of 
inflection  and  cause  a  roll-up  at  the  critical  layer  and  eventually  the  instability. 

The  influence  of  gas  density,  and  hence  aerodynamic  forces  on  the  interface,  has  been  widely  investigated 
experimentally.  Reitz  and  Bracco  [9]  observed  a  substantial  difference  in  the  atomization  mechanism  when 
the  liquid  jet  was  injected  in  different  gases  (i.e.  pNl  =  6  kg/m3  and  pXe  =  23  kg/m3).  Wu  et  al  [10] 
have  reported  a  change  in  droplet  size  for  primary  atomization  when  a  different  gas  density  was  tested  for 
the  same  liquid  jet.  Prom  this  basis  one  may  conclude  that  the  aerodynamic  interaction  at  the  surface  of 
the  jet  does  alter  wave  growth,  i.e.  that  aerodynamic  forces  are  of  sufficient  magnitude  to  contribute  to  the 
instability.  However,  Hoyt  and  Taylor[4]  did  not  see  differences  in  wave  structure  with  differing  external 
airflows  in  their  experiments.  The  axisymmetrically  disturbed  short  wavelength  observed  near  the  nozzle 
exit  is  present  regardless  of  magnitude  or  direction  of  the  air  velocity  as  noted  in  Fig.  l-(a).  Hoyt  and  Taylor 


59 


concluded  that  this  phenomenon  has  “no  discernible  effect  of  relative  air  velocity”.  In  spite  of  Hoyt  and 
Taylor’s  effort  in  this  regard,  their  issue  has  not  well  been  recognized  in  the  atomization  community. 

Shkadov’s  theory  states  that  the  velocity  profile  of  liquid  is  independent  of  gas  density,  which  seems 
contradictory  to  the  experimentally  observed  trend.  This  dilemma  is  probably  due  to  the  high  gas-liquid 
density  ratio  which  causes  the  growth  rate  due  to  the  aerodynamic  term  to  dominate  as  compared  to  the 
growth  rate  due  to  the  vortices,  thereby  invalidating  Shkadov’s  theory.  In  addition,  the  droplet  sizes  are 
strongly  dependent  on  secondary  atomization  at  high  jet  speed  and  high  gas  density  conditions.  For  many 
years,  the  notion  that  vorticity  at  the  nozzle  exit  being  responsible  for  the  atomization  has  been  over¬ 
shadowed  by  the  aerodynamic  linear  theory  and  the  experiment  which  showed  a  significant  difference  in 
spray  structure  with  the  change  in  gas  density.  Clearly  the  influence  of  the  gas  is  regime  dependent,  but  the 
fundamental  instability  mechanism  is  present  in  all  regimes. 

In  a  later  paper  by  Hoyt  and  Taylor  [11],  they  claimed  that  the  axisymmetrically  disturbed  wavelength 
shown  in  cavitating  flow  over  bluff  body  in  Fig.  2  is  similar  to  the  wavelength  seen  in  their  water  jet 
experiment.  Brennen  had  performed  linear  boundary  layer  instability  analysis  of  Rayleigh’s  equation  and 
provided  the  nondimensionalized  frequency,  7  =  0.175,  to  be  the  one  that  gives  the  maximum  amplification 
at  flow  separation  point.  Using  7  =  0.175  for  Hoyt  and  Taylor’s  case,  the  theoretically  predicted  wavelength 
was  A  =  (l/14.8)d  and  the  experimentally  observed  wavelength  was  A  =  (l/13.8)d  as  shown  in  Fig.  l-(b). 
The  comparison  between  theory  and  experiment  was  excellent. 

Based  on  the  evidence  that  Hoyt  and  Taylor  had  presented,  the  notion  that  boundary  layer  instability  is 
responsible  for  the  axisymmetrically  disturbed  waves  near  nozzle  exit  is  too  compelling  to  ignore.  For  this 
reason,  we  have  investigated  the  effect  of  boundary  layer  thickness  at  the  nozzle  exit  on  the  liquid  jet.  While 
the  entire  numerical  method  is  based  on  the  potential  theory,  it  may  sound  contradictory  to  model  the  effect 
of  viscosity  with  the  current  numerical  model.  However,  our  intention  is  to  flump’  the  vortices  at  the  nozzle 
exit  and  therefore  simulate  the  abrupt  change  in  boundary  condition  for  the  liquid  jet.  This  is  traditionally 
known  as  superposition  theory  and  is  somewhat  similar  to  what  has  been  known  as  the  Vortex-method’  (see 
Chorin  [12,  13,  14])  of  series  of  singular  vortices. 

A  complete  summary  of  the  boundary  layer  instability  is  discussed  and  contrasted  in  the  recent  study  by 
Yoon  and  Heister  [15].  Table  1  provides  the  summary  of  all  flow  regimes  mentioned  in  this  paper. 


Table  1:  Characteristics  of  various  flow  regimes 


Flow  Regime 

Rayleigh 

1st  WI 

2nd  WI 

Atomization 

Reitz’s  Criteria 

Weg  <  0.4 

0.4  <  Weg  <  12 

12  <  Weg  <  40 

Weg  >  40 

Reitz-Bracco’s  Recom.° 

Rayleigh 

Sterling-Sleicher0 

Kelvin-Helmholtz 

Reitz- Braccoc  (Taylor) 

Wavelength  Range*1 

A  >  4.12d 

0.73d  <  A  <  4.12d 

0.24d  <  A  <  0.73d 

A  <  0.73d 

Dominant  Forcee 

Capillary 

Cap.,  vise.,  dynm. 

Moderate  dynamic 

Strong  dynamic 

Dominant  Inst.  Mech7 

Capillary 

BLI 

BLI 

BLI  &  turbulence 

BLI  Effect 

None 

S  BL:  1 

Wiggles 

Finger-like  droplets,  lam.9 

M  BL:  Moderate  rollup  U  stretching 

Weak  BLI,  transitional 

L  BL:  Strong  rollup  &;  stretching 

turbulence  dominates 

WI:  wind  induced  “  applicable  eqn.  for  shear  layer  driven  instability 

b  for  ka  <  1.0  c  for  ka  -*•  oo 

d  based  on  Reitz-Bracco’s  recommendation  using  Reitz’s  criteria 


c,/  for  pi/pg  >  500;  shear  layer  driven  instability  is  negligible  [16] 

BLI:  boundary  layer  instability  S:  small  M:  medium  L:  large 

note:  the  range  of  S,M,L  cannot  be  determined  because  they  differ  at  different  nozzle  geometry  and  flow  condition 
9  this  flow  is  governed  by  inviscid  Orr-Sommerfeld  eqn. 


2  Modeling 

The  model  is  based  on  an  unsteady  axisymmetric  potential  flow  of  a  liquid  exiting  a  round  orifice  in  the 
absence  of  a  gas-phase  medium.  A  bound  ring  vortex  is  utilited  to  simulate  viscous  effects  associated  with 
vorticity  in  the  boundary  layer  formed  in  the  orifice  passage.  Carefully  controlled  experiments  have  shown 
a  reasonably  axisymmetric  structure  during  the  early  stages  of  the  free  surface  instability.  Fig.  3  provides 
a  schematic  representation  of  the  geometry  and  appropriate  nomenclature.  The  size  of  Rankine  vortex  [17] 
is  defined  as  Rc  which  will  be  discussed  in  detail  in  a  later  section.  A  vortex  ring  of  strength  r„  and 
overall  radius  f  is  assumed  to  lie  at  the  orifice  exit  plane.  A  computational  domain  represented  by  a  simple 
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cylindrical  column  of  length  zi  with  a  hemispherical  cap  is  selected  to  initialize  the  calculation.  Constant 
nodal  spacing,  As  is  employed  along  this  domain  and  nodes  are  added  as  the  jet  issues  forth  from  the  orifice. 
Fig.  4  compares  the  actual  flow  condition  and  the  current  superposition  modeling.  The  concentrated  vortices 
at  the  filament  vortex-ring  are  transported  to  the  free  surface  and  therefore  the  jet  surface  is  untable.  We 
choose  the  liquid  density,  p,  jet  average  exit  velocity,  U,  and  orifice  radius,  a  as  dimensions  in  the  problem. 

The  formulation  of  the  BEM  starts  with  the  integral  representation  of  Laplace’s  equation,  V20  =  0,  with 
(j>  being  the  velocity  potential.  Following  Liggett  and  Liu  [18],  the  integral  form  for  this  relation  is: 


dQ  =  0 


(i) 


where  <p{ft)  is  the  value  of  the  potential  at  a  point  f),  ft  is  the  boundary  of  the  3D  domain,  and  G  is  the 
free  space  Green’s  function  corresponding  to  Laplace’s  equation.  The  detailed  procedures  for  solving  Eq.  (1) 
using  BEM  is  available  in  the  literature  [19,  20,  21,  22], 

The  unsteady  Bernoulli  equation  provides  the  boundary  conditions  along  a  free  surface  [23].  This  condi¬ 
tion  provides  a  connection  between  the  inertial,  hydrostatic,  and  capillary  forces  at  the  interface.  Because 
the  surface  curvature,  k  depends  nonlinearly  on  the  surface  shape,  the  overall  expression  is  nonlinear.  Using 
the  nondimensionalization  described  previously,  the  appropriate  dimensionless  form  is: 
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Dt 


=  -\ut\2-ut-uv-  —  -Pg  + 


Bo 
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(2) 


where  Pg  is  the  dimensionless  gas-phase  pressure  (assumed  to  be  zero  in  the  present  studies),  and  We  and  Bo 
are  the  Weber  and  Bond  numbers  characterizing  the  flow:  We  =  pU2ajo  and  Bo  =  pga2/o.  This  Bernoulli 
Eq.  (2)  is  the  ‘nonlinear’  free  surface  boundary  condition  for  the  Laplace  equation.  This  equation  is  marched 
in  time  using  a  4th  order  Runge-Kutta  time  integration.  The  curvature  (k)  of  the  highly-distorted  surface  is 
determined  with  full  4th-order  accuracy  as  well  [24]  We  have  used  the  curvature  definition  from  Smirnov  [25], 
The  location  of  nodes  on  the  free  surface  (i.e.  z  and  r )  is  calculated  by  integrating  the  respective  velocity 
components  in  time. 

Contributions  from  the  ring  vortex  can  be  obtained  through  the  principle  of  superposition  for  potential 
flow.  Since  the  Laplacian  governing  equation  is  linear,  we  may  superpose  the  bulk  potential  flow  with  the 
potential  vortex-ring: 

4>t  =  <f>  +  <Pv  Ut  =  u  +  uv  (3) 

where  uv  is  the  induced  velocities  due  to  vortex  ring  that  can  be  obtained  from  the  Biot-Savart  law  [26,  27,  28]. 
The  entire  explicit  solution  of  the  filament  vortex-ring  has  recently  been  developed  [29],  Here  ( )t  represents 
the  general  or  ‘total’  solution  of  the  jet  flow.  The  solution  of  the  vortex  ring  can  be  obtained  by  direct 
evaluation  of  Biot-Savart  law  [28] .  The  entire  explicit  solution  of  the  filament  vortex-ring  has  recently  been 
developed  [29]. 

Nodes  are  repositioned  along  the  distorted  surface  using  cubic  splines  [30],  and  nodes  can  be  added 
(due  to  fluid  exiting  the  nozzle)  or  removed  (due  to  atomization  events)  without  user  intervention.  The 
current  ‘high-speed’  atomization  simulation  is  more  susceptible  to  numerical  instability  than  the  low-speed 
atomization  (i.e.  varicose  breakup  simulation  with  smooth  surface)  due  to  complex  surface  shape.  In 
addition,  the  ‘necking’  region  where  a  droplet  is  pinched-off  experiences  a  relatively  high  velocity  (i.e.  node 
velocities  reache  2  ~  7  times  that  of  the  jet  speed)  and  therefore  ‘node-crossing’  or  ‘not-simply-connected 
domain’  sometimes  occurs.  Thus  numerical  smoothing  is  crucial  to  prevent  the  numerical  instability  [31], 
We  have  chosen  the  filter  function  by  Longuet-Higgins  and  Cokelet  [31];  this  filter  is  designed  to  eliminate 
the  odd-even  mode  in  the  function  which  contains  the  highest  frequency,  known  as  the  Nyquist  Rate  [32]. 

The  centroid  of  the  vorticity  of  the  viscous  flow  is  regarded  as  the  center  of  the  vortex-ring.  The  definition 
of  the  centroid  of  the  vorticity,  weighted  in  the  radial  direction  is: 


C  n  Tudr 

Jr— 0 
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and  the  vorticity,  lj,  is  defined  as: 


where  u  and  v  are  the  velocities  in  axial  and  radial  direction  respectively.  Assuming  Qvjdz  »  0,  i.e.  parallel 
flow  at  the  orifice  exit  plane,  the  centroid  can  be  written  as: 


Iu(r=0)  TdU 

C:$d» 


(6) 


Substituting  the  definition  of  the  displacement  thickness,  Slt  [33]  into  Eq.  (6)  and  applying  integration  by 
parts  give  the  following  result: 

f  =  1  -  (7) 

can  be  approximated  using  a  Navier-Stokes  solution  of  the  internal  flow,  or  appropriate  analytical  methods 
such  as  Blasius  solution  [33]  for  a  flat-plate  or  Thwaites  [34]  equation  for  a  converging/diverging  nozzle. 
The  vortex  strength  is  defined  as  the  circulation  that  is  taken  about  any  path  enclosing  the  vortex-ring. 

Tv=Iu-dI  (8) 


where  u  is  the  internal  flow  velocity  of  the  injector  orifice  and  I  is  the  integration  path.  The  integration  is 
taken  in  counterclockwise  around  the  path.  Assuming  du/dz  =  0  and  a  no-slip  boundary  condition  at  the 
wall,  the  surface  integral  in  Eq.  (8)  can  be  written  as  follows  for  the  path  we  have  chosen: 

r„  =  Az  (9) 


where  A z  is  the  length  scale  which  is  comparable  to  the  most  dominant  wavelength,  A.  T„  is  always  positive 
and  will  induce  counterclockwise  motion  (this  is  based  on  the  upper  half  of  the  flow  going  from  left  to 
right).  Eqs.  (7)  and  (9)  uniquely  determine  the  location  and  strength  of  the  vortex  from  first  principles.  No 
additional  calibration  constants  are  used  in  the  formulation. 

It  is  possible  to  set  A z  =  A,  predicted  by  Brennen’s  result  [35]  since  the  circulation  causes  the  axisym- 
metrically  disturbed  wavelength  observed  at  the  nozzle  exit.  The  A  observed  in  Hoyt  and  Taylor’s  case  [11]  is 
a  function  of  the  momentum  thickness  S2,  scaled  by  the  parameter  7  =  0.175.  Presuming  a  high  contraction 
ratio  of  the  nozzle  reduces  the  turbulence  fluctuation,  Hoyt  and  Taylor  [11]  assumed  a  laminarized  flow  over 
a  flat  plat  and  therefore  they  utilized  the  Blasius  [33]  solution  to  approximate  the  momentum  thickness. 


Addition  of  the  Biot-Savart  Law  to  the  inviscid  jet  of  BEM  is  expected  to  cause  instability  at  the  free 
surface  that  eventually  forms  a  series  of  toroidal  ligaments  pinched  off  from  the  main  body  of  the  jet.  Using 
Gauss’s  divergence  theorem,  we  have  transformed  the  surface  integral  to  a  line  or  contour  integral  and 
therefore  have  obtained  the  cross-sectional  area  of  the  ligaments  as  well  as  the  centroids  of  the  area.  We 
have  obtained  the  volume  of  the  ligaments  using  the  theorem  of  Pappus-Guldinus  [36]  which  relates  a  volume 
of  revolution  to  its  generating  cross-sectional  area. 

While  the  current  model  is  based  on  an  axisymmetric  formulation,  actual  primary  atomization  is  a  3- 
dimensional  phenomenon  (see  Fig.  5-(a)).  The  model  result  of  pinch-off,  as  shown  in  Fig.  5-(b),  is  not  a 
droplet  but  a  vortex-ring  with  a  significant  amount  of  circulation  around  the  ring  surface.  The  circulation 
around  the  ring  surface  is  large  enough  to  cause  instability  in  circumferential  direction.  In  reality,  this  is 
the  secondary  instability  which  occurs  before  the  vortex-ring  pinch-off.  The  current  model  assumes  that 
droplets  are  formed  from  a  secondary  instability  on  annular  ligaments  shed  from  the  periphery  of  the  jet. 
This  amounts  to  a  decoupling  of  primary  and  secondary  instability  which  permits  the  axisymmetric  analysis 
of  the  jet  itself. 

Ponstein[l]  investigated  the  linear  stability  of  a  liquid  column  with  circulation  Tr  and  radius  ar.  Pon- 
stein’s  result  is  utilized  to  predict  the  dominant  wavelength  k  in  the  circumferential  direction: 
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where  Wer  -  pU2ar/cr.  Note  that  the  ring  radius,  ar,  is  the  non-dimensional  variable.  This  expression  is 
solved  to  determine  the  k  =  kmax  value  attributed  to  the  maximum  growth  rate,  u  for  a  given  ring  geometry 
and  circulation. 
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3  Results  and  Discussion 

3.1  Grid  Convergence  Study 

Hoyt  and  Taylor’s  case  is  used  for  the  grid  convergence  study  (i.e.  Welta  -  19057,  r*  =  0.99,  and  F*  =  0.139). 
Let  As  represent  the  grid  spacing  for  BEM  nodes.  While  Hilbing  [24]  mentioned  that  As  =  0.300  is 
fine  enough  to  resolve  the  low  speed  “Rayleigh’s  breakup”  where  waves  are  of  length  comparable  to  the 
orifice  diameter,  much  finer  grid  resolution  is  required  for  high  speed  atomization  where  the  wavelengths  are 
comparable  to  the  boundary  layer  thickness  at  the  orifice  exit.  For  this  reason,  the  grid  resolution  for  the 
present  studies  taxes  the  current  computational  capabilities  of  even  advanced  Linux-based  compute  clusters. 
In  Fig.  6,  it  is  shown  that  the  axial  location  for  the  first  ring  pinch-off  is  reasonably  insensitive  to  mesh  spacing 
for  As  <  0.030.  However,  grid  function  convergence  studies  indicate  that  a  smaller  mesh  spacing  is  required 
for  the  accurate  prediction  of  the  droplet  characteristics  in  the  atomization  regime.  About  3000  5000 

droplets  were  collected  for  each  run  for  statistically  reliable  data;  results  for  drop  statistics  are  shown  in 
Table  2.  The  Sauter  Mean  Diameter,  SMD,  (drop  whose  diameter  replicates  the  average  surface  area  of 
drops  in  the  population)  is  the  most  frequent  measure  used  in  the  atomization  field.  Table  2  shows  that  the 
Nd  per  ring,  its  Standard  Deviation,  and  SMD  are  converged  to  a  reasonable  accuracy  at  A  s/a  —  0.012. 
Thus  we  have  used  A s/a  =  0.016.  In  addition,  the  standard  deviation  of  N & /ring  does  not  change  much 
after  A  s/a  <  0.016.  Similarly,  the  time-averaged  droplet  velocities  (i.e.  uD  and  vD)  do  not  change  much 
either  for  this  mesh  spacing.  This  result  also  validates  the  pinch  criteria  employed  for  a  ligament  breakup:  a 
pinch-off  is  assumed  when  the  distance  between  binary  nodes  is  less  than  a  certain  tolerance,  e.  The  range 
of  the  tolerance  is  20  ~  70%  of  the  mesh  spacing.  The  ligament  size  is  also  insensitive  to  the  range  of  the 
pinch-off  criteria. 


Table  2:  Grid  convergence  test 


A  s/a 

N  o  /ring 

Stnd.  Dev. 

SMD/d 

*dIU 

vD/U 

UD/U 

0.050 

10.07 

4.68 

0.0995 

0.608 

0.442 

0.752 

0.040 

10.87 

8.71 

0.0962 

0.689 

0.430 

0.812 

0.030 

14.20 

11.33 

0.0841 

0.790 

0.479 

0.924 

0.020 

12.46 

8.46 

0.0729 

0.770 

0.419 

0.876 

0.016 

10.44 

6.65 

0.0623 

0.779 

0.404 

0.877 

0.012 

10.87 

6.44 

0.0588 

0.796 

0.411 

0.896 

It  is  uncertain  when  to  stop  the  simulation  since  the  jet  can  grow  indefinitely  depending  on  injection 
conditions.  For  the  simulations  conducted  to  date,  the  time  required  for  the  first  pinching  event  is  typically 
around  t*  &  1.7.  We  found  that  collecting  about  300  to  400  rings  provides  stastically  reliable  data.  This 
would  give  roughly  3000  ~  5000  droplets.  Thus  we  typically  stop  our  calculation  at  about  t*  &  5.0.  Table  3 
shows  that  the  solution  is  insensitive  when  t*  >  4.0. 


Table  3:  Effect  of  calculation  length  on  drop  statistics 


t* 

SMD  /d 

Nd 

N  D/ring 

Stnd.  Dev. 

2.0 

0.0628 

88 

7.97 

3.329 

3.0 

0.0623 

856 

10.44 

6.665 

4.0 

0.0635 

2582 

11.08 

7.936 

5.0 

0.0655 

5132 

11.48 

7.985 

3.2  Hoyt  and  Taylor’s  Case 

The  complete  simulation  of  the  Hoyt  and  Taylor’s  jet  [4]  is  shown  in  Fig.  7.  The  jet  structure  is  initially 
assumed  to  be  a  simple  cylinder  with  a  hemispherical  tip  as  shown  in  Fig.  3  and  its  evolution  is  simulated 
via  time  integration.  The  simulation  is  completed  at  t*  =  5.0.  A  slight  ‘swelling’  is  observed  at  t*  =  1.0 
and  a  fluctuation  of  the  jet  surface  is  seen  at  t*  >  2.0.  The  velocities  induced  by  the  bound  vortex  are  large 
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enough  to  penetrate  the  jet  surface  and  it  results  in  primary  atomization.  It  should  be  noted  that  most 
liquid  ligaments  pinching  from  the  jet  surface  are  in  the  ‘roll-up’  motion  in  counterclockwise  direction  while 
the  mean  velocity  of  the  ligament  is  in  the  stream-wise  direction.  Similar  structures  are  noted  in  Fig.  5- (a) 
in  a  closeup  view  of  the  Hoyt  and  Taylor  experiment.  The  mean  velocity  of  most  droplets  are  in  stream- 
wise  direction  as  droplets  motion  propagates  along  with  the  main  jet  stream,  the  most  dominant  convective 
source.  The  counterclockwise  roll-up  motion  is  a  strong  evidence  that  the  boundary  layer  instability  is  the 
fundamental  cause  of  the  primary  atomization.  The  counterclockwise  roll-up  motion  would  not  have  been 
observed  if  the  jet  were  turbulent  (see  Figure-(l)  of  Wu  et  al.  [16]).  The  mean  velocity  of  most  droplets  are 
in  stream-wise  direction  as  droplet  motion  propagates  along  with  the  main  jet  stream,  the  most  dominant 
convective  source. 

It  is  interesting  to  note  that  the  liquid  core  appears  naturally  as  a  consequence  of  the  calculation.  While 
the  current  model  is  based  on  the  axisymmetric  formulation  (2D),  the  real  instability  is  three-dimensional 
(3D).  It  is  obvious  that  the  loss  of  liquid  mass  of  the  model  prediction  is  noticeably  greater  than  that  of 
actual  3D  jet  and,  thus,  it  forms  the  liquid  core.  Another  possibility  for  accounting  for  the  difference  in  the 
actual  jet  and  the  model  jet  may  be  due  to  the  incomplete  simulation  of  the  actual  jet.  Presumably,  the 
entire  jet  should  be  streaming  down  far  away  from  near  nozzle  exit  region. 

All  droplets  pinched-off  within  t*  <  5.0  are  shown  in  Fig.  7.  It  is  seen  that  droplets  are  pinched-off  as 
soon  as  the  primary  instability  is  initiated;  this  may  seem  contradictory  while  experiment  shows,  in  Fig.  1- 
(a) ,  that  the  primary  instability  undergoes  the  transitional  process  to  result  in  the  primary  breakup  of  the 
secondary  instability.  The  difference  is  observed  because:  the  model  is  based  on  axisymmetrical  formulation 
and  does  not  include  viscous  effects.  During  the  transitional  region,  the  viscous  force  holds  the  fluid  together 
untill  the  vorticities  are  greater  to  win  over  the  viscous  force. 

Nevertheless,  the  model  predicts  the  droplet  diameter  (i.e.  d/ 10  <  D  <  dj 20)  which  is  well  within  the 
order  of  the  actual  SMD  of  Hoyt  and  Taylor’s  experiment  (i.e.  SMDexp  »  d/15.5),  as  shown  in  Fig.  13. 
Though  the  loss  of  mass  in  the  simulation  is  larger  than  that  observed  experimentally,  the  circulation  around 
the  pinch-off  ring  (rr)  is  large  enough  to  give  the  short  wavelength  which,  in  return,  accounts  for  the  less 
loss  of  mass  in  the  3D  actual  jet. 

The  pinch-off  droplet  velocity  vs.  its  location  is  plotted  in  Fig.  8.  It  is  shown  that  the  total  droplet 
velocity  (Uo  =  y/v%}  +  v2D)  increases  slightly  with  respect  to  the  pinch-off  axial-location  (z)  because  both 
ud  and  vo  increase  with  z.  This  is  consistent  with  Wu  et  al’s  observation  [16]  (see  Figure-(3)  of  Wu  et 
al.  [16]).  While  our  axial  velocity  (i.e.  uD  =  0.809)  is  in  excellent  agreement  with  Wu  et  al’s  result  (i.e. 
uD  =  0.760)  our  model  over-predicts  the  radial  velocity  (i.e.  vD  =  0.424)  as  compared  to  the  experimental 
result  (i.e.  vp  —  0.07).  The  difference  is  probably  due  to  (z)  the  recording  location  of  the  data:  while  our 
data  was  recorded  at  the  droplet  pinch-off  location,  Wu  et  al’s  data  was  recorded  at  a  few  nozzle  diameter 
away  from  the  centerline  (i.e.  r/d  ~  2  or  3)  where  droplet  had  already  experienced  drag.  This  could  possibly 
reduce  the  droplet  speed. 

Newton’s  2nd  law  is  applied  to  describe  the  motion  of  a  droplet  assuming  drag  to  be  the  only  external 
force  acting  on  a  droplet.  The  equation  of  motion  of  a  droplet  can  be  written  as: 
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Duo 

Dt 


—  2 Pg^oAp  |  Ug  Uo  I  ( ug  ud) 


(12) 


where  Ap  is  the  projected  area  of  a  droplet  ug  and  uo  are  gas  and  droplet  velocity,  respectively.  Co  is  the 
drag  coefficient  correlation  for  the  solid-sphere  given  by  Hwang  et  al  [37]  below: 


Co,s  = 


ReD  <  1000 
R eo  >  1000 


(13) 


where  ReD  =  UD/vair.  The  dynamics  of  a  droplet  of  Hoyt  and  Taylor’s  case  is  shown  in  Fig.  9.  It  indicates 
that  the  droplet  must  travel  a  distance  of  47  d  to  reduce  the  initial  radial  droplet  speed  (i.e.  vd  =  0.424  U )  to 
vd  =  0.2  U.  Thus  the  drag  alone  cannot  explain  the  difference  between  our  model  result  (i.e.  vd  =  0.424  U) 
and  the  experimental  observation  (i.e.  vp  =  0.07 17).  (ii)  It  is  possible  that  the  drag  coefficient  is  greater 
than  CD  -  0.518  for  a  solid-sphere  from  Hwang  et  al  [37]  because  CD  of  an  elliptically  distorted  droplet  in 
shape  is  usually  greater  than  that  of  a  solid-sphere.  (Hi)  Another  possibility  is  that  the  absence  of  viscosity 
in  our  model  may  have  induced  extremely  high  radial  velocity  during  the  pinch-off  event,  especially  at  the 
‘necking’  computational  node;  only  the  viscosity  effect  can  reduce  the  local  pinch-off  velocity. 
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We  have  observed  that  the  droplet  size  (D)  is  independent  (or  more  likely  random)  of  the  axial  location, 
a  smaller  D  is  observed  with  larger  radial  location.  We  have  enlarged  the  photo  of  Fig.  l-(a)  and  were  able 
to  confirm  that  the  droplet  closer  to  the  surface  is  larger  than  the  one  farther  away  from  the  surface. 

It  was  mentioned  in  the  previous  section  that  the  model  takes  advantage  of  Ponstein’s  Eq.  (11)  to  model 
the  instability  of  a  pinch-off  vortex-ring.  Using  Eq.  (11),  the  number  of  waves  or  droplets  {ND)  per  ring  is 
predicted  and  plotted  as  a  function  of  the  circulation,  Tr.  Their  relationship  is  parabolic  and  the  least  square 
fit  is  available  in  Fig.  10.  Thus  if  rr  is  known,  a  rough  estimate  of  ND  can  be  produced.  Fig.  11  illustrates 
the  relation  between  the  number  of  circumferential  waves  and  the  core  size  of  a  pinch-off  vortex-ring.  As  the 
size  increases,  a  larger  number  of  waves  appear.  This  is  exactly  the  opposite  phenomenon  of  the  elliptical 
instability  of  Widnall  [38]  and  Sullivan:  they  observed  a  larger  number  of  waves  when  the  core  size  decrease. 
It  should  be  noted  that  the  mechanism  of  surface  tension  driven  instability  (i.  e.  liquid  vortex-ring  in  air) 
is  different  from  the  shear  layer  driven  elliptical  instability  [39,  40,  41,  42]  (i.e.  liquid  into  liquid  or  gas  into 
gas). 

Hoyt  and  Taylor  [4,  11,  43]  reported  dramatic  changes  in  jet  configuration  even  with  small  (i.e.  10  ppm ) 
polymer  addition  to  the  water.  As  a  result  of  anisotropically  stretched  macromolecules  in  the  fluid  [44, 
45],  the  surface  drag  is  greatly  reduced  and  small-scale  surface  disturbances  are  eliminated.  However,  the 
polymer  addition  does  not  reduce  the  large-scale  disturbances.  In  fact,  Hoyt  et  al  [43]  reported  that  the 
polymer  solution  amplified  the  large-scale  disturbances  while  damping  out  the  small-scale  disturbances.  The 
instability  of  the  large-scale  waves  is  governed  by  Rayleigh’s  breakup,  A  =  4.51d,  in  sinuous  mode  even  for 
high  speed  jet.  Since  potential  superposition  theory  cannot  account  for  the  complex  governing  physics  of 
the  polymer  fluids,  only  non-polymer  jet  flow  is  simulated  in  this  report. 

3.3  Effect  of  Weber  Number  for  Fully  Developed  Flow 

We  may  consider  a  very  long  pipe  flow  with  constant  diameter  (i.e.  I/d  =  oo)  which  would  result  in  a 
fully-developed  flow  regardless  of  the  jet  speed,  U.  In  this  case,  the  jet  speed  determines  the  regime  of  the 
jet  instability  because  the  momentum  thickness  is  fixed  (and  therefore  the  circulation  amount  is  constant) 
for  all  speeds.  Thus,  more  atomization  is  expected  in  a  flow  with  the  high  jet  speed. 

Eventually,  the  boundary  layer  instability  disturbance  will  result  in  jet  breakup  regardless  of  the  jet  speed. 
If  U  is  small  that  Wei]a  =  100,  the  jet  similar  to  Rayleigh’s  breakup  is  observed  as  shown  in  Fig.  12  (note: 
no  atomization  is  observed  for  t *  <  5.0).  If  U  increases  up  to  Wel<a  =  1,000,  the  boundary  layer  instability 
effect  is  more  eminent;  ‘stretching’  is  observed  while  BLI  is  competing  with  a.  When  We^a  =  10, 000,  the 
surface  tension  force  is  greatly  reduced  and  therefore  the  atomization  event  occurs  at  a  faster  rate  (little 
competing  behavior  since  dynamics  force  is  dominant). 


Table  4:  Effect  of  Weber  number  on  drop  statistics.  Collected  Data  up  to  t*  =  5.0 


WeLa 

SMD/d 

No 

N  d /ring 

Stnd.  Dev. 

1,000 

0.1034 

1514 

3.03 

1.96 

10,000 

0.0733 

5172 

8.66 

5.87 

3.4  Effect  of  Jet  Speed 

Wu  et  al  [16]  provided  both  an  empirical  model  and  experimental  observations  for  the  Sauter  Mean  Diameter 
( SMD ).  Their  model  is  based  on  Kolmogorov  length  scale  [46]: 

SMD  133 

A  “  (14) 

Note  that  Wu  et  al  [16]  introduced  A  =  d/8  from  Hinze  [47].  The  above  expression  then  becomes: 

SMD  __  77.5  SMD  46.4 

d  ~  We°;J4  ~d~  ~  Wtfj 4  (15) 
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Here  SMD  is  defined  as: 


SMD  = 


Ellfi  Df 


(16) 


where  Nd  is  the  number  of  droplet  collected. 

It  is  well  known  that  the  droplet  size  varies  significantly  within  the  atomization  regime.  Wu  et  al  [16] 
reported  the  droplet  size  variation  with  U  for  turbulent  water  jet  into  air.  Hoyt  and  Taylor’s  experiment  had 
been  carried  out  for  the  Bernoulli  pressure  A P  <  60  psi;  no  result  with  higher  A P  is  reported  [4,  11,  48]. 
However,  we  hypothesized  the  increase  in  A P  for  Hoyt  and  Taylor  jet.  The  result  has  taken  the  jet  speed 
up  to  U  =  40  m/s  which  corresponds  to  A P  «  116  psi  or  slightly  higher  for  accounting  some  pressure  loss 
within  the  nozzle.  The  final  result  for  Hoyt  and  Taylor’s  case  [4]  is  shown  in  Fig.  13.  Using  the  methodology 
employed  in  the  prior  section  with  no  calibration  constants,  the  model  predicts  the  Sauter  Mean  Diameter, 
SMD,  with  reasonable  accuracy.  As  shown  in  Fig.  13,  there  is  a  steep  gradient  at  jet  speed  around  U  «  20  m/s 
for  Wu  et  al  [16]  turbulent  jet  experiment.  Our  model  result  overlaps  with  that  obtained  by  Hoyt  and  Taylor. 
It  is  interesting  to  observe  that  Wu  et  al’s  data  is  also  similar  to  our  result  and  that  of  Hoyt  and  Taylor  at 
about  U  «  20  m/s.  While  Wu  et  al  noted  this  as  the  region  of  ‘uncertainty’,  it  is  possible  that  the  roll-up 
motion  was  competing  with  the  turbulence  and  thus  the  perceptible  effect  of  the  roll-up  motion  appears  as 
shown  in  Fig.  13. 

As  U  >  20  m/s,  differences  between  the  calculations  and  the  experiments  merge.  It  is  known  that  linear 
analysis  [49]  overpredicts  the  droplet  size  (by  less  than  20%)  because  it  neglects  the  satellite  droplet  mass 
due  to  the  nonlinear  effect  [50,  51]  which  yields  multiple  crests  per  wavelength.  Ponstein’s  Eq.  (11)  is  a  linear 
analysis  and  thus  it  also  may  tend  to  overpredict  the  droplet  size.  However,  20%  or  smaller  SMD  difference 
does  not  explain  the  difference  we  see  in  Fig.  13  at  higher  jet  speed.  It  is  due  to  the  fundamental  difference 
between  the  boundary  layer  instability  jet  and  the  turbulent  jet:  Boundary  layer  instability  jet  is  scaled 
by  the  momentum  thickness  [35]  and  the  turbulent  jet  is  scaled  by  the  Kolmogrov  length  scale,  lk  or/and 
turbulence  eddy  characteristics  length  of  kinetic  energy,  U  [16].  Wu  et  al  derived  the  empirical  formula  using 
the  ‘surface  kinetic  energy’  argument  which  gives  SMD  scaled  by  ~  l/U1Aa.  Thus  the  governing  length  scale 
(i.e.  lk  and  h)  decrease  significantly  at  about  U  ~  20  m/s.  On  the  other  hand,  SMD  of  the  boundary  layer 
instability  jet  is  scaled  by  S2:  SMD  ~  1/U0-5  and  thus  its  change  with  respect  to  U  is  relatively  moderate 

as  shown  in  Fig.  13.  . 

It  is  noted  that,  in  reality,  cavitation  might  occur  at  the  sharp  corner  of  the  nozzle  entry  during  this 
high  speed  flow  ( U  ~  40  m/s)  because  it  is  known  that  the  cavitation  magnitude  increases  as  the  jet 
speed  increases  [52].  In  addition,  the  jet  may  become  turbulent  if  Re&2  >  Res2tCru  due  to  increase  in  U 
(i.e.  Res2  oc  U1/2).  In  Fig.  14,  the  jet  becomes  more  stable  (or  less  atomization)  as  U  increases.  This  is 
consistent  with  Rupe’s  observation  [5].  As  U  increases,  the  nozzle  exit  velocity  profile  becomes  flatter  which 
contains  less  circulation.  One  may  imagine  the  limiting  case  of  such  condition:  U  — ►  oo  ( Red  -*  oo)  and 
thus  inviscid  flow  with  the  perfectly  flat  velocity  profile  which  contains  no  circulation  at  the  nozzle  exit.  In 
this  limiting  case,  the  jet  is  unconditionally  stable  unless  it  is  perturbed  by  other  instability  mechanisms.  If 
the  circulation  at  the  nozzle  exit  increases,  the  jet  is  more  unstable.  However,  the  greater  stability  of  the 
higher  speed  BLI  jet  structure  does  not  necessarily  mean  that  the  higher  speed  jet  produces  lesser  number 
of  droplets.  In  fact,  the  circulation  of  the  pinch-off  ring  for  the  higher  speed  jet  is  greater  and  therefore  it 
results  in  the  smaller  droplet  size  and  more  total  number  of  droplets  up  to  a  certain  limit  (note:  the  limit 
has  not  yet  been  determined  in  this  study).  In  Fig.  15,  it  is  shown  all  droplet  statistics  (i.e.  total  number  of 
droplets,  Nd,  and  its  value  per  ring,  Nd/ ring,  and  its  standard  deviation)  increase  with  increasing  U .  This 
indicates  that  the  droplet  statistics  becomes  more  unstable  as  U  increases. 

In  summary, 

•  The  SMD  of  the  boundary  layer  jet  is  scaled  by  1/U 0  5  while  that  of  the  turbulent  jet  is  scaled  by 
1/U148. 

•  The  smaller  circulation  with  increasing  U  at  the  nozzle  exit  yields  more  stable  jet  structure. 


•  On  contrary,  the  circulation  around  the  pinch-off  ring  contains  the  larger  circulation  with  increasing 
U.  Thus  it  yields  the  smaller  droplet  size  and  more  unstable  droplet  statistics. 
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3.5  Effect  of  Exit  Plane  Boundary  Layer  Thickness 

Changing  the  ratio  of  nozzle  length-to-diameter  ( l/d )  is  essentially  the  same  as  changing  the  boundary  layer 
thickness  for  the  constant  diameter  pipe  flow.  The  effect  of  l/d  for  the  jet  breakup  has  been  investigated  by 
many  researchers  [53,  54,  55,  9,  56]. 

For  the  Rayleigh  and  1st  wind-induced  regime  (Weg  <  2.55),  Sterling-Sleicher’s  experimental  data  shows 
that  the  breakup  (or  jet)  length,  L,  decreases  with  increasing  nozzle  length,  l.  This  indicates  that  the  velocity 
profile  with  thicker  boundary  layer  breaks  up  faster  due  to  larger  roll-up  motion  at  the  critical  layer.  They 
had  also  observed  the  important  effect  of  the  velocity  profile  relaxation  as  they  reduced  the  aerodynamic 
effect  by  decreasing  the  ambient  pressure  from  0.98  to  0.2  atm. 

For  the  2nd  wind-induced  regime,  McCarthy  and  Molloy  [56]  also  investigated  the  effect  of  l/d  on  at¬ 
omization  mechanism  as  shown  in  Fig.  16.  The  jet  is  shown  up  to  zjd  <  70  and  the  jet  speed  remains 
constant.  The  flow  operating  condition  was:  Liquid:  60%  glycerol  and  40%  water  by  weight,  pi  =  103 , 
Pi  =  11  cP,  ai  =  0.0669^,  U  =  20  d  =  2.54mm,  Weiia  =  781,  Red  =  4750,  and  Weg  =  18.7  for  l/d=0™l, 
5,  and  10.  We  can  see  that  the  flow  is  in  the  low  part  of  the  2nd  wind-induced  regime  and  thus  no  definite 
dominating  instability  mechanism  exist  (see  Table  1).  The  boundary  layer  instability  is  competing  with 
other  instabilities  such  as  capillary  and  viscous  force.  As  to  which  force  should  be  dominant  is  determined 
by  the  boundary  layer  thickness  during  this  flow  regime.  For  instance:  the  instability  was  not  observed  for 
l/d  =  0.  For  l/d  —  1,  the  relaxation  length  was  about  lr/d  «  4.7.  The  axisymmetrically  disturbed  waves 
(due  to  primary  instability)  were  seen  at  this  point  but  the  viscous  force  is  large  enough  to  pull  the  surface 
structure  intact.  Only  wiggles  appear  on  the  surface  due  to  the  competition  between  vorticity  and  viscosity. 
For  l/d  =  5,  the  relaxation  length  is  increased  to  lr/d  ~  5.7  and  the  primary  atomization  is  observed  (due 
to  secondary  instability).  For  l/d  =  10,  the  relaxation  length  is,  again,  increased  to  lr/d  «  6.4.  SMD  does 
not  seem  to  change  much  from  the  previous  case  while  the  atomizing  motion  became  more  amplified  due  to 
larger  roll-up  motion.  The  increase  in  relaxation  length  with  larger  boundary  layer  appears  to  be  reasonable; 
it  takes  longer  distance  for  the  critical  layer  to  develop  the  roll-up  motion  for  the  thicker  boundary  layer! 
The  roll-up  motion  causes  the  jet  to  be  atomized  when  the  circulation  is  large  enough  to  win  against  other 
competing  forces  such  as  viscous  or/and  capillary  force.  In  Fig.  17- (top),  the  ‘stretching’  is  observed  due 
to  capillary  force  when  the  jet  is  atomized.  In  Fig.  17-(bottom),  the  computational  result  for  the  case  is 
shown  where  the  similar  stretching’  of  the  capillary  force  is  observed.  For  the  high  speed  jet  like  that  of 
Hoyt  and  Taylor  [4,  11],  the  effect  of  viscosity  and  capillary  force  is  of  little  importance  because  the  jet  is 
nearly  inviscid.  In  fact,  the  large  scale  motion  of  Hoyt  and  Taylor’s  jet  is  governed  by  the  Rayleigh’s  inviscid 
analysis  [49]  which  concludes  the  most  dominant  wavelength  to  be  A  =  4.5 Id. 

If  the  viscosity  is  strong  enough  to  damp  out  the  ‘thin’  or  small  boundary  layer  instability  and  thus 
competing  behavior,  wiggles  on  the  surface  are  observed.  It  is  known  that  the  viscosity  generally  damps 
out  the  instability  (even  though  sometimes  the  instabilities  arise  due  to  the  viscosity,  such  as  Tollmien- 
Schlichting  waves  [7]).  In  other  words,  the  circulation  of  the  ‘thin’  boundary  layer  is  not  strong  enough  to 
result  in  the  atomization  in  the  2nd  wind-induced  regime  (though  ‘thin’  boundary  layer  in  the  atomization 
regime  is  strong  enough  to  result  in  the  atomization  due  to  the  absence  of  v  and  small  <x).  The  model  cannot 
predict  this  behavior  due  to  the  absence  of  viscosity,  as  noted  in  Table  1.  The  model,  which  assumes  the 
inviscid  flow,  is  more  unstable  than  the  actual  flow  due  to  the  absence  of  viscosity.  Perhaps  the  inviscid  flow 
with  rotationality  is  intrinsically  unstable  unless  the  surface  tension  is  relatively  large  enough  to  suppress  the 
rotationality  or/and  dynamic  force  of  the  fluid.  However,  the  model  predicts  that  ‘stretching’  atomization 
occurrs  with  the  ‘thick’  boundary  layer  within  the  1st  or/and  2nd  wind-induced  regime,  as  shown  in  Fig.  17- 
(bottom).  In  this  case,  the  atomization  occurs  in  a  ‘stretching’  pattern  while  the  roll-up  motion  is  enhanced 
by  the  surface  tension  force,  which  tends  to  keep  the  surface  stable.  McCarthy  and  Molloy’s  experiment  is 
summarized  in  Table  5,  which  shows  an  excellent  agreement  between  the  wavelength  predicted  by  Brennen’s 
boundary  layer  analysis  and  the  experimental  results. 

For  the  atomization  regime  {Weg  >  40),  Reitz  and  Bracco  [9]  show,  in  Fig.  6  in  their  paper,  that  the 
cone  angle  increases  with  increasing  l/d.  This  trend  is  evident  only  in  the  high  gas  density  environment  and 
with  high  cavitation  in  the  internal  flow.  It  is  well  known  that  increasing  gas  density  results  in  higher  cone 
angle.  It  is  probably  that  the  radial  velocity  of  droplets  ( v& )  increases  with  the  gas  density,  which  makes  the 
cone  angle  greater;  Wu  et  al  [10]  had  observed  that  v&  increases  with  increasing  gas  density.  This  suggests 
that  increasing  l/d  for  the  Reitz  and  Bracco’s  case  might  have  had  the  similar  effect  (i.e.  causing  the  higher 
vo)  to  cause  the  higher  cone  angle. 
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Table  5:  Summary  for  McCarthy  and  Molloy’s  experiment  [56] 


WKM 

81/d 

T  [Hz] 

Ac 

A  exp 

1 

mwnmm 

BjQtgsl 

d/2.8 

5 

wm&m 

d/1.4 

mm 

47477 

1/12.66 

1/32.82 

145 

7197 

1.1  d 

1.0  d 

a  where  Rex  =  Ul/v  b'c  Brennen’s  Eqs.  (??)  and  (??),  respectively 
note:  Blasius  solution  is  used  for  <Si  and  S2  estimation 


In  another  experiment  for  the  atomization  regime,  Hiroyasu  et  al  [54]  (i.e.  water  jet  into  air)  shows 
that  the  increase  in  l/d  causes  the  breakup  length  L  to  be  smaller  (This  is  consistent  with  Sterling  and 
Sleicher’s  observation  even  though  the  flow  regime  is  different).  However,  beyond  the  maximum  value  at 
about  U  «  60  m/s,  their  behavior  shows  no  clear  trend.  In  addition,  this  trend  could  not  be  observed  when 
Hiroyasu  et  al  [54,  55]  changed  the  ambient  pressure  from  0.1  to  3.0  MPa;  aerodynamic  effect  is  greater. 

We  have  learned  that  the  influence  of  l/d  is  difficult  to  investigate  according  to  these  researchers  work. 
The  reason  is  that  there  are  many  other  mechanisms  (i.e.  cavitation,  turbulence,  and  aerodynamics)  that 
may  be  coupled  with  the  effect  of  l/d.  In  addition,  the  influence  of  l/d  becomes  less  perceptible  when  the 
flow  regime  approaches  to  the  atomization  regime.  However,  we  have  noticed  one  clear  trend  (even  though 
we  may  draw  different  conclusion  from  the  work  of  these  researchers):  “Thicker  boundary  layers  are  more 
unstable” . 

For  this  reason,  we  have  investigated  the  effect  of  boundary  layer  thickness:  Three  cases  for  different 
boundary  layer  thickness  are  conducted  using  the  current  model  (i.e.  l/d  =  1,  2,  and  00).  The  configuration 
of  the  test  cases  for  the  Hoyt  and  Taylor’s  nozzle  is  shown  in  Fig.  18.  Only  the  boundary  layer  thickness 
changes  while  everything  else  remains  constant  (i.e.  Reynolds  and  Weber  numbers).  It  is  shown  in  Fig.  19 
that  the  increase  in  S2  results  in  a  larger  fluctuation  at  the  early  stage  of  the  atomization  (even  though 
its  effect  fades  away  as  t*  increases  in  Fig.  22).  A  good  explanation  for  the  faster  instability  growth  with 
larger  S2  is  depicted  in  Fig.  20:  as  the  critical-layer  moves  away  from  the  parallel  line,  it  allows  a  larger 
roll-up  motion  which  results  in  larger  fluctuation  of  the  jet  structure.  Note  that  Lin  [57,  58]  shows  that  the 
critical-layer  is  located  at  the  point-of-inflection  of  the  velocity  profile.  Table  6  summarizes  the  results  of 
three  cases  considered.  A  slight  increase  in  SMD  is  seen  with  increasing  £2-  This  is  due  to  larger  droplets 
arising  from  larger  roll-up  motion.  The  larger  fluctuation  gives  more  dispersive  No  distribution  as  its  RMS 
value  increases  from  Case-1  to  Case-3.  This  dispersive  behavior  with  larger  S2  is  seen  in  Fig.  21.  One 
may  suspect  the  limiting  case  of  such  behavior  with  the  fully  developed  flow:  S2  is  a  maximum  under  such 
conditions.  A  fully  developed  flow  is  likely  to  contain  substantial  turbulence  fluctuation  as  the  laminarized 
flow  is  limited  only  in  the  entry  region.  Turbulent  flow  is  expected  to  have  the  most  dispersive  behavior  in 
the  droplet  distribution.  As  time  progresses,  it  is  difficult  to  observe  the  difference  between  Cases  1  to  3 
as  shown  in  Fig.  22.  In  addition,  little  change  in  SMD  is  found  (see  Table  6).  This  is  the  similar  difficulty 
that  the  previous  researchers  [54,  55]  had  experienced.  Even  though  our  current  model  is  not  influenced  by 
other  effects  (i.e.  turbulence,  cavitation,  or/and  aerodynamics),  we  have  to  agree  that  the  influence  of  the 
boundary  layer  itself  is  difficult  to  observe  when  the  spray  is  in  the  atomization  regime.  However,  it  is  clear, 
as  shown  in  Fig.  23,  that  Case  3  of  'thicker’  boundary  layer  results  in  more  unstable  condition  with  more 
atomizations. 


Table  6:  Model  predictions  for  different  boundary  layer.  Collected  data  up  to  t*  =  5.0 


Case 

l/d 

81/d 

82/d 

SMD/d 

nd 

Up 

vd 

N  0 /ring 

Stnd.  Dev. 

1 

1.0 

1/200 

1/518 

1/15.27 

3884 

0.809 

0.424 

11.48 

7.98 

2 

2.0 

1/141 

1/367 

1/15.19 

5879 

0.817 

0.428 

11.59 

8.25 

3 

00 

1/6 

1/24 

1/14.81 

8133 

0.811 

0.457 

11.69 

8.40 
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4  Conclusions 


^  fully  nonlinear  model  lias  been  developed  to  simulate  primary  atomization  caused  by  boundary  layer 
instability  using  superposition  of  a  ring  vortex  with  a  potential  jet  flow.  The  axisymmetric  model  employs 
a  boundary  element  methodology  in  which  the  vorticity  in  the  boundary  layer  at  the  orifice  exit  is  used  to 
determine  ring  vortex  strength  and  radial  location  at  the  orifice  exit  plane.  Annular  ligaments  are  pinched 
off  the  surface  in  this  case;  a  secondary  linear  instability  analysis  due  to  Ponstein  is  used  to  predict  the 
fractionization  of  the  ligaments  into  individual  droplets. 

The  SMD  of  the  model  result  agrees  well  with  the  actual  droplet  size  of  Hoyt  and  Taylor’s  experiment. 
The  result  of  the  current  model  is  also  compared  with  the  experimental  data  of  Wu  et  al.  [16],  The  comparison 
confirms  that  the  current  model  predicts  the  droplet  size  satisfactorily. 
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Figure  2:  Typical  cavitating  flow  over  bluff  body  (ogive).  Experimental  image  by  Brennen  [35].  Printed 
under  the  Permission  of  Journal  of  Fluid  Mechanics. 


Figure  3:  Schematic  of  the  initial  jet  geometry  indicating  computational  nodes  and  the  axisymmetric  ring 
vortex  at  the  orifice  exit  plane. 
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Figure  4:  Comparison  of  the  actual  physical  flow  condition  with  the  superposition  model  in  simulating 
boundary  layer  relaxation  downstream  of  the  orifice  exit  plane. 


(a) 


Figure  5:  (a)  Closeup  of  the  actual  Hoyt  and  Taylor’s  water  jet  [4].  Printed  under  the  permission  of  Journal 
of  Fluid  Mechanics,  (b)  The  closeup  of  the  model  result  for  Hoyt  and  Taylor’s  water  jet. 
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Figure  8:  Model  result:  velocity  contour  in  axial  and  radial  direction  for  the  Hoyt-Taylor  case  [11] 


80 


Figure  9:  Dynamics  of  a  droplet  of  Hoyt  and  Taylor’s  case:  d  =  6.35  mm,  D  =  d/15.5,  U  =  21  m/s, 
ReD  =  UD/uair  =  589,  uD  =  0.809  U,  vD  =  0.424  U,  pi  =  999  kg/m3,  pg  =  1.23  kg/m3,  and  CD  =  0.518.  ’ 
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Figure  10:  Prediction  of  cicumferential  wave  number  (or  number  of  droplet)  due  to  circulation  around  the 
rotating  ring  pinched-off  from  the  main  liquid  stream. 


Figure  11:  Number  of  waves  vs.  core  thickness  of  pinch-off  vortex-ring 
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Figure  12:  Effect  of  Weber  number  on  atomization  of  fully-developed  flows. 
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Figure  13:  Sauter  mean  diameter  comparison  for  U  ~  20m/, 


Figure  14:  Effect  of  jet  speed  on  jet  surface  structure  of  Hoyt-Taylor’s  jet.  Due  to  larger  circulation  contained 
in  the  velocity  profile,  more  unstable  jet  structure  appears  in  the  lower  speed  case. 
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Figure  15:  Effect  of  jet  speed,  £7,  on  drop  statistics 
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(a)  Experiment 


(b)  Model  result 


Figure  17:  Liquid:  60%  glycerol  and  40%  water  by  weight,  p  =  103^,  p  =  llcP,  a  =  0.0669^£,  U  =  20—, 
d  =  2.54mm,  We  =  781,  Pe  -  4750.  Comparison  between  experiment  [56]  and  model  results.  The  black 
circles  represent  the  location  and  the  relative  size  of  the  pinehed-off  droplets.  Printed  under  the  permission 
of  Elsevier  Science. 
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Figure  18:  Pressure  distribution  of  Hoyt-Taylor’s  nozzle  [11]  for  three  different  nozzle  lengths.  Note  that 
‘Case-3’  is  the  fully-developed  flow  (i.e.  I/d  =  oo). 
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Roll-up  motion  becomes 
larger  as  82  increases 


Velocity  profile  Location  variation  of 

point-of-inflection  (critical-layer) 


Figure  20:  Schematic  of  rollup  motion  at  critical  layer  for  the  three  cases. 
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t  =3.00 


t  =4.00 


Figure  22:  Boundary  layer  effect  fades  away  as  t*  increases 
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8  Appendix  D  -  Coaxial  Injector  Simulations 


Kim*,  B  and  Heister,  S.  D.,  “Two-phase  Modeling  of  Hydrodynamic  Instabilities  in  Coax¬ 
ial  Injectors”. 
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Abstract 

The  hydrodynamic  instability  of  coaxial  injector  with  recessed  inner  fluid  posts  have  been  investi¬ 
gated  using  a  homogeneous  flow  method  with  pseudo-density  model  in  a  time-dependent,  viscous 
calculation.  The  present  study  focuses  on  unsteady  self-oscillation  mode  in  the  coaxial  injector. 
The  Kelvin-Helmholtz  instability  mechanism  due  to  velocity  discontinuity  between  gas  and  liquid 
phase  is  investigated  as  a  source  of  unsteadiness  which  could  contribute  to  combustion  instabili¬ 
ties.  A  2-D  analogue  is  investigated  in  these  initial  studies.  A  series  of  parametric  studies  have 
been  completed  to  assess  the  amplitude  and  frequency  of  oscillations  under  a  variety  of  design 
conditions.  Effect  of  the  thickness  of  the  inner  post  has  been  studied  as  well. 


Introduction 

In  many  coaxial  injectors,  liquid  flowing  down  a  central  post  is  atomized  by  a  high-velocity  gas 
passing  around  the  outer  annulus.  In  many  applications,  the  liquid  post  is  submerged  somewhat 
from  the  orifice  exit  plane  to  provide  flame  holding  in  combustion  systems  such  as  liquid  rocket 
engines.  High-frequency  oscillation  of  the  jet  spray  has  been  observed  by  various  researchers’ 
experiments  of  coaxial  injectors  with  the  submerged  liquid  post.  1.2,3,6—8,21  xhese  disturbances 
could  couple  to  the  dynamics  of  the  jet  breakup  process  and  potentially  provide  amplification 
of  oscillations  within  the  combustion  chamber.  Bazarov  1>2  dubbed  this  phenomenon  as  ‘Self- 
Oscillation’  and  suggested  this  as  a  cause  of  decreased  combustion  efficiency  and  a  source  of  high 
amplitude  noise  during  combustion.  Combustion  instabilities  of  this  nature  can  have  severe  impact 
on  the  performance  of  the  engine  and  can  in  some  cases  lead  to  catastrophic  failures. 

Hutt  and  Rocker 4  investigated  the  high  frequency  combustion  instability  associated  with  coax¬ 
ial  injectors.  They  classified  the  instability  phenomena  in  the  chamber  as  injection-coupled  and 
intrinsic  mechanism.  The  injection  coupling  implies  chamber  pressure/temperature  variation  as  a 
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key  contributor  in  the  change  of  flow  dynamics  through  the  injector.  In  the  other  hand,  the  intrinsic 
mechanism  occurs  in  the  flowfield  due  to  its  own  flow  dynamics  with  negligible  feed  system  effect. 
However,  it  should  be  noted  that  injection  coupling  is  never  independent  of  the  intrinsic  subpro¬ 
cesses,  such  as  atomization,  propellant  heatup,  vaporization,  and  mixing  because  these  processes 
determine  the  relationship  between  the  injector  response  and  the  chamber  response. 

Bazarov  x’ 2  studied  the  self-oscillation  phenomena  along  with  the  self-pulsation  mode  of  jet 
instability  in  coaxial  injector.  He  postulated  that  the  self-oscillation  occurs  when  the  gas-liquid 
interaction  forms  a  cavity  inside  nozzle,  leading  to  jet  swirling  around  the  nozzle  exit.  The  self¬ 
pulsation  of  the  liquid  jet  mixed  with  the  gas  flow  depends  on  the  pressure  drop  at  liquid  and  gas 
phase,  correlating  with  the  time  of  liquid  propagation  through  the  injector  nozzle  5.  In  addition  to 
external  disturbances  from  combustion  chamber  and  feed  system  during  engine  operation,  injectors 
can  generate  self-pulsation  under  certain  conditions  by  its  own  intrinsic  unsteadiness  leading  to 
random  modeification  of  the  spray  formation  process. 

Mayer  and  his  research  group  have  done  a  significant  amount  of  work  on  the  coaxial  injector 
in  terms  of  the  combustion  instability.  Mayer  and  Krulle  6  investigated  coaxial  flow  mixing  phe¬ 
nomena  in  terms  of  chamber  pressure  variation,  density/velocity  ratio  changes  and  surface  tension 
effect.  By  increasing  chamber  pressure  gas  density  is  increased,  magnifying  the  aerodynamic  in¬ 
teraction  between  the  liquid  and  gas  phases  and  resulting  in  a  faster,  finer  atomization.  Increasing 
surrounding  gas  velocity  also  leads  to  an  increase  of  surface  wave  growth  and  to  macroscopic  in¬ 
stabilities  of  the  liquid  jet.  They  claimed  the  initiation  of  the  jet  surface  deformation  was  due  to 
internal  liquid  turbulence  delivering  energy  transform  in  forms  of  eddy  structures,  approximately 
a  size  of  10-30%  of  the  LOX  post  diameter.  In  other  works  of  Mayer  7,s,  the  coaxial  injector  flow 
was  studied  under  firing  engine  condition  at  supercritical  chamber  pressure  (higher  than  5  MPa). 
The  study  revealed  a  remarkable  difference  between  subcritical  spray  formation  and  the  supercrit¬ 
ical  injection  and  mixing.  At  subcritical  condition,  the  liquid  jet  shows  similar  flow  pattern  to  the 
cold  flow  test  forming  ligaments  off  the  liquid  jet  surface  and  producing  droplets  before  evapora¬ 
tion.  Upon  approaching  and  exceeding  supercritical  pressure,  droplets  no  more  exist  and  the  liquid 
jet  rapidly  dissolves.  The  flame  from  combustion  chamber  was  anchored  at  the  tip  of  LOX  post 
by  flow  recirculation  eddies  serving  as  flame  holder  for  steady-state  combustion.  The  asymmet¬ 
ric  flow  oscillation  was  also  reported  in  all  experiments,  but  the  source  of  the  oscillation  was  not 
clearly  stated. 

Instability  mechanisms  in  coaxial  injectors  were  also  investigated  experimentally  under  non¬ 
combusting  conditions  by  Glogowski  et  al.  9,10  Their  experimental  results  showed  that  for  the 
coaxial  injector  with  a  liquid  oxygen  (LOX)  post  recessed  into  the  fuel  annulus,  the  injector  tran¬ 
sitioned  into  a  condition  of  resonance  characterized  by  a  whistling  noise.  Significant  modification 
to  the  overall  structure  of  the  spray  due  to  the  strong  acoustic  coupling  between  injector  hydrody¬ 
namics  and  spray  formation  was  also  found.  Without  the  recessed  region,  the  injector  operation 
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produced  a  near  resonance  condition  with  a  lower  amplitude  whistling  noise  but  did  not  make  con¬ 
siderable  change  in  spray  structure.  The  recirculation  by  reverse  flow  near  the  recessed  LOX  post 
exit  was  also  showed  in  this  study.  The  effort  to  measure  flow  properties  within  the  recessed  region 
was  not  successful  due  to  optical  constraints  and  the  small  spatial  and  temporal  scales  involved. 

Eroglu  and  Chigier 11  focused  on  the  wave  characteristics  of  liquid  jet  from  coaxial  air-blast 
injector.  They  measured  frequency  and  wavelength  of  the  jet  issuing  from  the  injector,  and  found 
two  dominant  wave  type;  spanwise(dilational)  and  stream-wise(sinuous)  waves.  The  spanwise 
wave  showed  at  a  low  relative  jet  velocity  between  gas  and  liquid-phase  and  the  stream-wise  wave 
showed  at  a  high  relative  jet  velocity.  Average  wavelengths  decreased  with  liquid  and  gas  velocity. 
The  frequency  band  of  the  jet  oscillation  increased  with  the  liquid  jet  velocity. 

Mansour  and  Chigier  12  also  conducted  similar  research  on  the  liquid  sheet  instability  issuing 
from  the  two-dimensional  air-assisted  nozzle.  The  results  showed  the  same  pattern  with  the  Eroglu 
and  Chigier’s  study.  However,  Mansour  and  Chigier’s  study  probed  higher  velocity  cases,  and 
confirmed  that  the  frequency  of  the  liquid  sheet  oscillation  increased  with  coflowing  gas  velocity. 
This  also  indicates  that  the  aerodynamic  interaction  between  gas  and  liquid  flow  is  the  dominant 
factor  for  the  flow  oscillation. 

In  general,  for  the  flow  at  high  velocity,  most  of  researchers  agree  that  the  principal  source 
introducing  instability  to  the  jet  is  from  aerodynamic  forces  arising  from  the  interaction  of  the 
liquid  jet  with  the  surrounding  gas  flow.  Reynolds  and  Weber  numbers  are  generally  very  high 
in  these  atomizers  and  aerodynamic  forces  are  several  orders  of  magnitude  larger  than  capillary 
forces.  The  interaction  between  the  liquid  and  gas  phases  mainly  comes  from  different  velocities 
of  each  phase.  The  velocity  discontinuity  in  a  homogeneous  fluid  results  wave  growth  on  the 
interface,  which  is  a  common  Kelvin-Helmholtz  instability. 

Most  of  the  previous  research  works  summarized  here  mainly  contributed  to  the  investiga¬ 
tion  of  flow  phenomena  outside  of  injector  after  the  injection.  None  of  them  addressed  the  inner 
flow  structure  in  a  recessed  region  of  the  coaxial  injector,  especially  using  numerical  simulation. 
Since  the  upstream  flow  in  a  injector  provides  the  initial  condition  for  the  entire  spray  atomization 
process  and  combustion  chamber  acoustic  characteristics,  research  work  on  this  area  is  highly  de¬ 
sirable.  In  this  paper,  the  2-D  analogue  (a  liquid  sheet  in  a  channel)  of  the  coaxial  injector  will 
be  studied  numerically  to  characterize  amplitude  and  frequency  of  oscillations  of  the  inner  liquid 
sheet.  The  model  is  described  briefly  in  the  following  section;  results  of  parametric  studies  are 
then  summarized  in  the  latter  part  of  the  paper  along  with  the  result  of  finite  LOX  post  thickness 
modeling. 
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Linear  Analysis 


Liquid  sheet  instabilities  have  been  analyzed  by  many  researchers  over  the  past  several  decades. 
Among  those  investigations,  Squire’s  17  analysis  was  among  the  first  fundamental  studies  on  the 
instability  of  inviscid  liquid  film.  Based  on  the  Squire’s  analysis,  Dombroski  and  Johns  18  inves¬ 
tigated  the  aerodynamic  instability  of  viscous  liquid  film  subject  to  symmetric  and  asymmetric 
waves.  Li  and  Tankin  19  also  studied  viscous  liquid  sheets  and  provided  an  analysis  for  aerody¬ 
namic  and  viscosity-enhanced  instability  modes.  However,  most  previous  efforts  have  focused  on 
a  moving  liquid  sheet  in  a  semi-infinite  stagnant  gas.  Since  the  jet  inside  coaxial  injector  encoun¬ 
ters  higher  velocity  gas  phase  and  the  outside  nozzle  wall  provides  a  finite  region  for  the  gas,  the 
mathematical  formulation  for  this  special  case  needs  to  be  derived.  In  this  section,  an  analytical 
approach  based  on  the  linear  theory  is  presented  for  the  coaxial  injector  case. 

It  should  be  noted  that  the  approach  based  on  the  linear  stability  theory  can  only  evaluate  the 
beginnings  on  the  instability.  As  the  amplitude  of  the  disturbance  grows,  the  linear  assumptions 
soon  become  invalid.  In  addition,  the  presence  of  a  finite-length  domain  makes  the  connection  to 
an  infinite  sheet  stability  analysis  somewhat  tenuous:  i.e.,  the  injection  exit  plane  is  not  treated  in 
periodic  boundary  conditions  typically  used  in  linear  studies.  Even  though  the  theoretical  approach 
doesn’t  allow  us  to  describe  the  jet  instability  in  a  complete  structure,  it  provides  a  good  starting 
point  for  predicting  unstable,  wave-type  behavior  of  the  jet  with  disturbances. 

Consider  a  two-dimensional  liquid  sheet  with  density  pi,  surface  tension  a  and  uniform  thick¬ 
ness  2 h,  moving  at  velocity  Ui  through  an  inviscid  moving  gas  medium  of  density  pg  at  velocity 
Ug  as  illustrated  in  Figure  1 .  The  coordinates  are  chosen  such  that  the  direction  of  the  z-axis  is 
parallel  to  the  direction  of  the  velocity  Ui  and  Ug,  and  the  y- axis  is  normal  to  the  liquid  sheet  with 
its  origin  located  at  the  mid-plane  of  the  liquid  sheet.  Let  the  relative  velocity  A U  represents  the 
velocity  difference  between  gas  and  liquid  phase,  Ug  -  Uh  that  is,  the  reference  frame  is  attached 
to  the  liquid  phase.  Finally,  let  h  represent  the  sheet  thickness  and  nh  represent  the  distance  to  the 
channel  wall,  where  n  measure  the  number  of  sheet  thicknesses  to  the  wall. 

For  anti-symmetrical  disturbances,  the  displacements  of  corresponding  points  on  the  two  sur¬ 
faces  are  equal  in  magnitude  and  in  the  same  direction.  Hence,  the  two  interfaces  are  regarded  to 
have  the  following  form  : 


y  =  ±h  +  i)  (1) 

r]  =  r]oe^t+ikx)  (2) 

where  y  =  ±h  are  the  equilibrium  positions  of  the  two  interfaces,  i.e.  the  position  with  no  distur¬ 
bances;  rjo  is  the  initial  amplitude  of  the  wave,  and  is  taken  to  be  much  smaller  than  the  half-width 
of  the  sheet,  h;  k  is  the  wave  number  of  the  disturbance,  and  k  =  2ir/\,  where  A  is  the  wavelength 
of  the  disturbance;  to  =  ujr  +  ioj,  is  a  complex  variable.  The  real  part  ^represents  the  rate  of 
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growth  or  decay  of  the  disturbance;  its  imaginary  part  w*  is  2ir  times  the  disturbance  frequency; 
and  —oji/k  is  the  wave  propagation  speed  of  the  disturbance. 

Let  <f>  =  (f)g  +  A Ux  be  the  velocity  potential  of  the  gas  phase  where  the  <j)g  is  the  perturbation 
potential  and  fa  is  the  perturbation  velocity  potential  of  the  liquid  phase.  We  will  neglect  viscosity 
for  this  analysis  since  the  Reynolds  numbers  are  quite  high  in  practical  devices.  Therefore,  the 
Laplace  equation  governs  the  process: 


d2<j>  d2(j) 

dx2  dy2 


(3) 


where  <f)  is  the  velocity  potential  with  fa  and  (f>g  for  the  individual  liquid  and  gas  phases,  respec¬ 
tively.  The  linearized  kinematic  boundary  condition  for  liquid  phase  is 


dfa  _  dr] 
dy  dt 


(4) 


and  for  gas  phase 


=  +  atj^L 

dy  dt  dx 


(5) 


which  are  to  be  satisfied  at  y  —  ±h.  The  dynamic  boundary  condition  from  unsteady  Bernoulli’s 
equation  gives  the  following  relation  for  the  pressure  at  the  interface: 


Pi 


dx  , 


Since  the  pressure  induced  by  surface  tension  a  is. 


(6) 


Pi 


(7) 


to  2nd  order  in  tj.  Here,  R  is  the  radius  of  curvature  of  the  interfaces.  Considering  the  disturbances 
given  in  Equation  (2),  it  is  assumed  that  fa  and  <f>g  take  the  following  form: 


fa  =  F(y)e^t+ik^  (8) 

<t>9  =  G{y)e^t+ikx^  (9) 

Equation  (8)  and  (9)  can  be  solved  by  performing  a  Fourier  analysis.  The  velocity  potential  for  the 
liquid  phase  becomes 

fa  =  T]0  tanh(kh)e^t+lkx )  (10) 

and  the  velocity  potential  for  the  gas  phase  is: 


1 fig  — 


r/o  cosh(ky)  —  tanh(nkh)sinh(ky) 
k  sinh(ky)  —  tanh(nkh)cosh(ky)  J 


(w  +  ikAU)e^t+ikx) 


(11) 
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Substitution  of  Equation  (10),  (1 1)  and  (7)  into  Equation  (6)  for  y  =  h  leads  to  the  following 
dispersion  relation  between  the  complex  growth  rate  uj  and  the  disturbance  wave  number  k  : 

u?{pitanh{kh)Dh  —  pgNh )  —  u(2ikAUpgNh )  +  k2  AU2  pgNh  +  okzDh  =  0  (12) 

where 

Ni j  =  cosh(kh)  —  tanh(nkh)  sinh(kh)  (13) 

Dh  =  sinh(kh )  —  tanh(nkh)cosh(kh )  (14) 

Equation  (12)  represents  a  complex  quadratic  equation  for  growth  rate  a;  as  a  function  of  wave 
number  k. 

The  dispersion  relation  can  be  nondimensionalized  using  liquid  density  pi,  liquid  phase  velocity 
Ui  and  the  jet  thickness  h  as  dimensions.  Assuming  Weber  number  We  »  1  due  to  the  high  liquid- 
phase  velocity,  the  normalized  dispersion  relation  becomes: 

u2(tanh(k)Dh  —  eN /,)  —  uj{2ikAUtNh)  +  k2eAU2Nh  =  0 

and  since 

Nh  cosh{k)  -  tanh(nk)sinh(k )  _  1  -  tanfe(nA:)ton/i(fc) 

Dh  ~  sinh(k)  —  tanh(nk)cosh(k )  tanh(k)  —  tanh(nk ) 

by  letting  Nh/Dh  be  ©,  the  dispersion  relation  becomes  : 

uj2(tanh{k)  —  e@)  —  co(2ikAUeQ)  +  k2eAU2Q  =  0 

Since  the  condition  for  instability  is  where  the  real  part  of  growth  rate  is  positive,  solving  this 
dispersion  relation  for  growth  rate  and  evaluating  the  effect  of  changing  the  three  main  parameters 
(h,  e,  A U)  can  provide  us  basic  idea  of  how  the  liquid  sheet  will  behave. 

In  order  to  verify  this  analytic  solution,  infinite  distance  between  the  liquid  sheet  surface  and 
the  nozzle  wall  is  applied,  i.e.,  n  oo.  In  this  limit,  the  hyperbolic  tangent  term,  tanh(nk)  has 
an  approximate  value  of  1,  and  the  Equation© 6)  becomes: 

g*  =0==  l~t*nh(k)  _  (18) 

Dh  tanh(k )  —  1 

Under  the  long  wave  assumption.  Equation  (17)  becomes  : 

w2(l  +  e)  +  u)(2ikAUe)  -  k2eAU2  =  0  (19) 

This  result  is  identical  to  Squire’s  result  in  the  long  wave  limit.  The  parameter  ©  really  measures 
the  effect  of  the  finite  channel  height.  Therefore,  solving  Equation  (17)  identifies  the  influence  of 
the  finite  channel  height  of  interest  in  this  work. 


(15) 

(16) 

(17) 
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Figure  2  shows  the  nondimensional  growth  rate,  ur  for  anti-symmetrical  disturbances  when 
the  relative  velocity,  AU  varies.  Here  the  density  ratio,  e  is  0.01  and  the  channel  width  to  sheet 
thickness  ratio,  n  is  5.  Cases  of  Ug/Ui  =  2,4  and  8  are  investigated.  As  shown  in  the  result, 
the  liquid  sheet  instability  grows  rapidly  as  the  gas  velocity  grows.  Due  to  the  high  speed  of  the 
gas  phase  surrounding  liquid  sheet,  the  Kelvin-Helmholtz  instability  governs  the  behavior  of  the 
liquid  sheet.  While  the  growth  rate  increases  indefinitely  as  the  wave  number  grows,  viscosity  and 
surface  tension  in  actual  flows  will  reduce  the  growth  rate  in  large  wave  number  region. 

Figure  3  illustrates  the  effect  of  gas/liquid  density  ratio  on  growth  rate.  Here  AU  =  3  and 
n  =  2  for  all  cases,  and  density  ratio  e  values  of  0.1,  0.01  and  0.001  are  considered.  As  gas  phase 
density  increases,  the  gas  imparts  more  momentum  to  the  liquid  surface.  Therefore,  the  growth 
rate  of  disturbance  increases  when  the  density  ratio  gets  higher. 

The  effect  of  channel  width  to  jet  thickness  ratio,  n  is  shown  in  Figure  4.  The  density  ratio,  e 
is  0.01  and  the  relative  velocity,  AU  is  3,  and  n  values  of  2,  5,  and  10  are  considered.  This  result 
shows  that  the  instability  grows  rapidly  as  the  sheet  gets  thicker  (channel  gets  smaller).  However, 
in  the  small  wave  number  region  (k  <  0.07)  the  growth  rate  exhibits  the  opposite  pattern.  The 
presence  of  the  cavity  wall  has  a  stabilizing  effect  on  the  very  long  wave  instabilities.  This  behav¬ 
ior  is  most  important  in  assessing  the  growth  rate  response  since  the  long  wavelength  instabilities 
are  presumed  to  be  the  most  harmful  to  overall  injector  performance  in  combustion  system.  There¬ 
fore,  one  can  conclude  that  the  presence  of  wall  can  either  be  a  stabilizing  or  destabilizing  effect 
depending  on  the  wavelength  of  the  instability. 


Modeling  Description 

A  two-dimensional  incompressible,  unsteady,  viscous  flow  solver  has  been  developed  utilizing  a 
finite  volume  implementation  of  the  Marker  and  Cell  discretization  method.  The  current  model  is 
based  on  homogeneous  flow  two-phase  treatment  in  which  the  single  phase  Navier-Stokes  equa¬ 
tions  are  solved  using  a  fictitious  “pseudo”  density  which  varies  in  amplitude  between  the  liquid 
and  gas  extremes.  This  provides  a  mechanism  to  compute  the  local  droplet  number  density  with¬ 
out  having  to  solve  the  flow-field  around  all  the  individual  droplets.  The  single  fluid  model  can 
be  achieved  by  assuming  locally  homogeneous  flow  (LHF)  in  which  the  relative  velocity  and  tem¬ 
perature  between  two-phases  are  small  enough  in  comparison  to  variation  of  the  overall  flow  field 
that  is  to  be  predicted.  Under  the  LHF  assumption,  and  by  providing  a  proper  constitutive  relation 
for  the  pseudo-density  of  the  homogeneous  flow,  the  model  is  able  to  handle  the  two-phase  flow 
with  less  computational  resources  than  traditional  two-fluid  modeling.  However,  capillary  forces 
are  not  resolved  using  this  approach  because  the  interface  is  not  known  as  part  of  the  solution 
methodology.  For  the  high  Reynolds  and  Weber  numbers  characterizing  the  injectors  of  interest, 
this  simplification  is  deemed  appropriate.  The  development  of  the  homogeneous  fluid  model  is 
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discussed  in  detail  in  prior  works  13  15 . 

The  two-dimensional,  viscous,  unsteady,  Navier  Stokes  equations  are  expressed  in  the  follow¬ 
ing  form: 


—  + 

dpv 

=  0 

dt 

dx 

dy 

dpu 

dpu2 

dpuv 

dP 

_  d 

du 

d 

du 

dt 

+  dx 

+  dy 

dx 

dx 

dy 

^  dy 

dpv 

dpuv 

dpv2 

dP 

_  d 

dv 

d 

dv 

dt 

+  dx 

+  dy 

+  dy 

dx 

f‘^  + 

dy 

^ dy 

The  Lagrangian  form  of  the  continuity  equation  is  also  required  : 


Dp  .du  dv . 


(20) 

(21) 


(22) 


(23) 


Because  of  the  two-phase  treatment,  the  viscosity  can  vary  spatially.  According  to  Kubota  et 
al.  16 ,  the  viscosity  of  mixture  can  be  written: 


p  =  apg  +  (1  -  a)pi 


(24) 


where  pg  and  pi  are  the  gas  and  liquid  viscosities,  and  a  is  the  void  fraction.  Since  the  non- 
dimensional  pseudo-density  is  volume  fraction  of  mass  per  unit  cell  volume,  the  Equation  (24)  can 
be  written  as: 

p(p)  =  ppi  +  (1  -  p)pg  (25) 


This  mixture  viscosity  is  substituted  back  into  Equation  (21)  and  (22)  for  non-dimensionalization. 
The  channel  width,  liquid  inflow  velocity  and  liquid  density  are  chosen  as  dimensions  in  nondi- 
mensionalizing  the  equations.  Rearranging  the  equations  for  flux  calculation  through  cell  faces 
yield  the  following  momentum  equations: 
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(27) 


where 


1  P  ,  Pg0--P)  1 


_  ^ -  (28) 

Re*  Rei  pi  Reg 

which  includes  the  both  liquid  and  gas  phase  viscous  effect  in  one  term.  The  second  and  third  terms 
in  left  hand  side  of  the  Equation  (26)  and  (27)  are  the  main  terms  for  calculation  of  momentum 
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fluxes  in  x  and  y  directions  in  the  code.  Here  Ret  and  Reg  represent  Reynolds  numbers  of  liquid 
and  gas-phase  respectively. 

Since  p  is  a  non-physical  variable,  an  additional  constitutive  relation  is  required  in  place  of 
an  equation  of  state  which  would  normally  close  the  set  of  governing  equations.  In  this  case, 
we  envision  a  group  of  droplets  convecting  through  a  gas  media.  If  we  assume  all  droplets  are 
of  the  same  size  and  neglect  fission  or  fusion  processes,  then  the  local  number  density  sets  the 
pseudo-density  for  a  given  computational  cell.  Consequently,  one  can  conclude  from  the  above 
assumptions  that  the  density  evolves  in  time  by  simply  considering  a  Lagrangian  tracking  of  the 
droplet  field  as  specified  by  a  continuity  equation  as: 


Dp 

Dt 


dp  dp  dp 

di  +  uai  +  %  =  0 


(29) 


This  equation  is  basically  a  statement  that  the  droplet  mass  is  invariant  along  a  path  line  in  the 
fluid.  By  taking  account  for  the  flow  direction,  the  pseudo-density  is  updated  based  on  mass  flux 
calculation  as  follows : 


Pn+1  =  Pn  ~  ^  (mout)  -  fmin) 

Vc  [  V  /  x,y  \  /  x,y_ 


(30) 


where  Vc  is  volume  of  a  given  cell  and  m  is  mass  flowrate  corresponding  x  and  y  directions.  The 
flux  of  the  mass  through  a  given  cell  are  calculated  by  a  standard  upwind  scheme.  While  this 
provides  a  locally  first  order  solution  in  the  region  adjacent  to  the  discontinuity  in  density,  it  pro¬ 
vides  for  a  stable  approach  to  account  for  the  large  density  gradients  near  the  interface.  Since  the 
interface  region  will  undoubtedly  contain  droplets  in  these  very  high  convective  environments,  a 
shock-type  density  discontinuity  would  not  be  consistent  with  a  homogeneous  flow  representation. 
The  scheme  is  second-order  accurate  with  the  exception  of  the  points  adjacent  to  the  interface,  and 
convergence  is  verified  in  the  next  section.  Figure  5  illustrates  the  basic  geometry  of  a  coaxial  in¬ 
jector.  The  computational  domain  is  the  recessed  region  inside  the  nozzle,  and  the  structured  mesh 
employed  100  x  200  grid  points  in  transverse  (y)  and  axial  (x)  directions  respectively.  Exponential 
stretching  is  applied  in  the  transverse  direction  to  enhance  resolution  near  the  walls. 

For  boundary  conditions,  the  liquid  and  gas-phase  velocities  are  defined  at  inflow  boundary. 
The  pressure  is  extrapolated  with  zero  gradient  for  inflow  boundary.  No  slip  conditions  for  veloc¬ 
ities  are  defined  on  both  walls.  Finally,  constant  pressure  condition  is  set  for  outflow  boundary 
while  the  velocities  are  extrapolated  at  this  location. 

A  series  of  grid  convergence  studies  have  verified  that  this  mesh  (100  x  200)  is  adequate  to 
resolve  the  unsteady  sheet  flow.  When  the  velocities  for  the  liquid  sheet  and  surrounding  gas  flow 
in  upper  and  lower  sides  are  balanced  in  exact  symmetry,  the  flow  field  reaches  a  steady  state  after  a 
certain  period  of  time.  This  steady  state  condition  provides  us  an  opportunity  for  relatively  simple 
grid  function  convergence  test. 


105 


A  center  point  in  the  computational  domain  has  been  picked  for  a  flow-representative-position. 
Time  histories  of  density,  pressure  and  liquid  velocity  have  been  observed  until  the  flow  solution 
reaches  to  the  steady  state.  Figure  6  shows  the  each  flow  property  history  at  the  center  point  for 
three  grid  sizes:  100  x  50, 200  x  100  and  400  x  200  in  x  and  y  axis  respectively.  Result  shows  that 
the  200  x  100  grid  provides  accuracy  comparable  to  the  fine  grid.  In  Figure  7  and  8,  the  density 
and  velocity  distribution  over  y  axis  at  the  nozzle  exit  are  depicted  after  reaching  steady  state. 
The  coarse  grid  reaches  steady  state  when  t  =  3.7,  and  the  medium  and  fine  grids  achieve  steady 
state  at  t  =  5.15,  t  =  4.75  respectively.  These  figures  also  show  the  good  agreement  between  the 
medium  grid  and  fine  grid  while  the  coarse  grid  still  has  some  errors  as  compared  to  the  other  two 
grids.  In  overall,  the  200  x  100  grid  is  expected  to  generate  flow  solution  in  a  good  resolution  while 
reducing  running  time  in  a  significant  order.  A  typical  calculation  takes  three  to  four  hours  for  the 
200  x  100  grid  on  Pentium  II  450MHz  machine.  The  400  x  200  grid  runs  about  48  hours  even  on 
64-bits  RISC  Alpha  chip  (500MHz)  machine  until  the  flow  solution  reaches  to  the  steady  state. 

Since  it  takes  substantially  long  time  until  the  liquid  sheet  develops  the  instability,  an  artificial 
disturbance  is  introduced  in  order  to  initiate  the  oscillation  at  ealier  stage.  This  has  been  done  via 
setting  up  unbalanced  velocities  on  upper  and  lower  side  of  liquid  surface  grid  points  at  initial  time 
step(f  =  0.0)  only.  The  result  shows  that  a  certain  level  of  small  disturbance  can  cause  the  violent 
oscillation  of  the  liquid  sheet  inside  the  channel  after  a  short  transient  time.  Since  the  flow  field 
is  solved  by  a  system  of  elliptic-parabolic  governing  equations,  the  numerical  disturbance  causing 
the  oscillation  is  initial  condition  dependent.  However,  by  no  means  of  measuring  the  magnitude 
of  the  artificial  disturbance  quantitatively  at  the  time  the  oscillation  is  initiated,  a  direct  comparison 
between  the  cases  with  and  without  the  initial  disturbance  is  not  possible.  The  best  one  can  do  at 
the  moment  is  to  make  sure  that  the  code  predicts,  at  least  qualitatively,  all  characteristics  of  the 
flow. 


Numerical  Results 

A  baseline  case  is  presented  to  provide  the  reader  with  insight  into  the  unsteady  jet  oscillation 
inside  a  coaxial  injector  nozzle.  The  injector  schematic  is  shown  in  Figure  5.  The  channel  width, 
Dq,  is  assumed  to  be  0.005m  and  the  liquid  sheet  width,  Dif  is  0.001m.  The  length  of  recessed 
region  L  is  0.01m  which  is  twice  the  channel  diameter.  The  thickness  of  liquid  injector  structure 
is  neglected  since  the  computational  domain  is  outside  of  the  liquid  injector  and  the  gas  and  liquid 
flows  interact  each  other  right  after  coming  out  of  the  exit.  The  gas  and  liquid  phase  velocities 
are  80m/s  and  20 m/s  respectively.  The  Reynolds  numbers  are  Reg  =  2.22  x  105  for  gas  and 
Rei  —  9.98  x  104  for  liquid.  Here,  density  of  water  has  been  used  for  the  liquid  phase  and  a 
gas/liquid  density  ratio  of  0.01  is  assumed  to  be  realistic  with  high  pressure  injector  operation. 

The  overall  behavior  of  the  instability  is  shown  in  Figure  9  which  depicts  pseudo-density  con- 
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tours  at  various  times  during  the  calculation.  Here,  the  outermost  contour  is  for  pg  =  0.01  and  the 
innermost  contour  is  for  pi  =  1.0.  After  some  chaotic  oscillations  during  early  phases  of  the  calcu¬ 
lation,  the  sheet  enters  a  quasi -periodic  oscillation.  The  velocity  streamlines  and  velocity  profiles 
at  the  exit  are  shown  in  Figure  10. 

By  plotting  the  location  of  maximum  density  at  the  exit  plane  as  a  function  of  time,  the  char¬ 
acter  of  the  oscillation  becomes  apparent  as  seen  in  Figure  1 1  which  shows  a  very  high  amplitude 
oscillation  (channel  walls  are  at  r  =  0,  1)  as  evidenced  in  the  pseudo-density  contours  in  Figure  9. 

The  oscillation  is  also  quasi-periodic;  by  taking  the  Fast  Fourier  Transform  of  the  signal  in 
Figure  1 1,  we  can  identify  its  fundamental  frequencies.  Results  of  this  calculation  are  shown  in 
Figure  12.  The  peak  near  /  =  0  is  attributed  to  the  transient  during  the  start  of  the  calculation. 
The  first  peak  with  frequency  of  0.6836  indicates  the  fundamental  frequency,  and  the  second  peak 
corresponds  to  the  first  harmonic  mode  at  frequency  of  1.3867.  The  development  of  multiple  har¬ 
monics  takes  longer  than  for  a  single  tone.  There  is  some  activity  at  higher  harmonics,  but  at  a 
much  lower  energy  level  than  the  primary  tone  and  the  first  harmonic.  Under  the  nondimensional- 
ization  employed,  a  frequency  of  2.0  would  correspond  to  the  time  it  takes  a  liquid  fluid  element 
to  traverse  the  channel  according  to  the  definition  of  Strouhal  number.  The  primary  tone  would 
correspond  to  a  period  of  roughly  three  times  the  channel  transit  time. 

Converting  these  nondimensional  frequency  values  into  physical  ones,  the  primary  harmonic 
is  at  2740  Hz  and  the  second  harmonic  is  near  5600  Hz.  These  are  frequencies  within  the  range 
of  acoustic  modes  within  liquid  rocket  engine  combustion  chambers.  A  series  of  experimental 
result  by  Eroglu  and  Chigier 11  also  showed  that  the  frequencies  of  jet  oscillation  of  airblast  coaxial 
injector  at  similar  water-air  velocities  fall  into  the  same  frequency  band.  In  principal,  the  numerical 
analysis  does  support  the  conclusion  that  jet  instability  in  the  submerged  region  could  reinforce 
instabilities  in  the  combustion  chamber.  However,  the  jet  submergence  is  much  larger  here  than 
that  typically  used  in  rocket  engine  injectors.  For  this  reason,  a  parametric  study  was  initiated  to 
classify  the  instability  over  a  range  of  design  and  operating  conditions. 

Parametric  Study 

A  parametric  study  has  been  conducted  in  order  to  evaluate  the  influence  of  gas/liquid  density  ratio, 
velocity  ratio,  sheet  thickness,  sheet  submergence,  and  Reynolds  number.  Results  are  compared 
by  plotting  the  amplitude  of  the  oscillation  (basically  the  difference  between  the  upper  peaks  and 
lower  peaks  in  Figure  1 1)  as  a  function  of  time  for  the  various  parameters  investigated. 

Density  ratio  effects 

Four  simulations  at  various  gas/liquid  density  ratios  with  all  other  inputs  fixed  (Ug /Ui  =  4,  L/D  = 
2,  h  =  0.2,  Re  =  105)  have  been  conducted  in  order  to  assess  the  influence  of  this  parameter.  The 
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oscillation  amplitude  histories  for  these  cases  are  summarized  in  Figure  13  for  gas/liquid  density 
ratios,  e  =  0.005,  0.01,  0.02,  and  0.1.  The  primary  frequency,  /,  of  each  case  is  also  noted  in  the 
legend  in  Figure  13.  Results  indicate  that  both  the  amplitude  and  the  frequency  of  the  oscillation 
increase  with  gas  density.  Physically,  this  is  attributed  to  the  fact  that  the  gas  is  capable  of  imparting 
more  momentum  to  the  liquid  when  it  is  of  higher  density.  This  behavior  is  completely  consistent 
with  the  linear  Kelvin-Helmholtz  theory. 

Effect  of  gas/liquid  velocity  ratio 

Three  different  velocity  ratio  cases  have  been  studied  for  other  conditions  consistent  with  the  base¬ 
line  case  described  previously.  Here,  Ug/Ui  =  2,  4,  and  8  were  investigated;  results  are  depicted 
in  Figure  14.  The  case  with  Ug/Ui  =  2  shows  very  little  activity  as  the  sheet  deviates  little  from  its 
initial  height.  The  case  with  Ug/Ui  -  8  setting  shows  a  chaotic  oscillation  which  essentially  en¬ 
compasses  the  entire  width  of  the  channel.  In  this  case,  it  appears  that  the  sheet  breaks  up  because 
of  the  unusually  high  gas  velocity  and  the  liquid  fragments  disintegrated  from  the  liquid  core  keep 
hitting  walls  and  generate  non-periodic,  high  frequency  vibration  rather  than  wave  type  oscillation. 
For  this  reason,  no  frequency  value  is  reported  for  this  condition. 

Physically,  the  gas  dynamic  pressure  which  drives  the  instability  scales  as  the  square  of  the  gas 
velocity.  The  results  obtained  in  this  study  is  consistent  with  this  general  nonlinear  behavior. 

Effect  of  varying  sheet  thickness 

Geometric  variations  were  also  investigated  in  the  studies.  Figure  15  depicts  the  influence  of 
the  thickness  of  the  liquid  sheet  at  other  conditions  corresponding  to  the  baseline  case.  Sheet 
thicknesses  ( h )  corresponding  to  20,  40,  60,  and  80%  of  the  channel  width  were  considered  in  the 
study.  Physically,  this  would  correspond  to  thicknesses  of  1,  2,  3,  and  4  mm  in  the  baseline  injector 
design.  In  general,  the  oscillation  frequencies  were  not  strongly  dependent  on  the  thickness  of  the 
sheet.  The  amplitude  of  the  oscillations  tends  to  decrease  with  increasing  h,  presumably  due  to  the 
fact  that  the  liquid  inertia  grows  with  the  sheet  thickness. 

Linear  analyses  indicate  that  the  presence  of  the  walls  has  a  destabilizing  effect  which  is  the 
opposite  of  the  trends  noted  in  the  nonlinear  calculations.  The  long  wave  instabilities  which  govern 
the  calculations  are  indeed  damped  by  the  presence  of  channel  walls.  High  frequency,  short  wave 
instabilities  are  not  resolved  in  the  present  model  since  capillary  forces  are  neglected  and  the 
pseudo-fluid  treatment  does  not  support  a  sharp  discontinuity  at  the  interface.  This  area  would 
require  further  study  with  a  separate  flow  two-phase  model  in  order  to  more  fully  resolve  this 
issue. 
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Channel  length  effects 

The  effect  of  liquid  injector  submergence  is  addressed  in  Figure  16.  Here,  the  liquid  injector  was 
recessed  at  four  separate  distances,  1,  1.5,  2,  and  3  channel  heights  (L/D  =  1,  1.5,  2,  3)  with 
other  conditions  remaining  identical  to  the  baseline  configuration.  When  the  channel  width  is  the 
same  as  the  recess  length  (L/D  =  1),  there  is  little  activity,  and  a  straight  sheet  comes  out  of 
the  exit.  As  the  submergence  length  is  increased,  the  L/D  =  1.5  case  shows  the  sheet  starts  to 
oscillating  but  the  magnitude  of  the  oscillation  is  smaller  than  the  baseline  case  and  damps  out 
gradually  as  time  goes  on.  At  L/D  =  2,  large  amplitude  oscillations  are  noted. 

At  L/D  =  3,  the  oscillation  becomes  unpredictable  with  varying  magnitude.  In  this  case,  the 
liquid  sheet  breaks  up  before  it  reaches  to  the  exit,  and  the  liquid  fragments  from  the  disintegration 
near  the  exit  are  shot  out  at  somewhat  irregular  intervals.  The  radical  change  of  the  amplitude 
of  this  case  can  be  explained  by  this  behavior.  In  general,  these  trends  are  in  agreement  with  the 
linear  theory  since  increased  submergence  increases  the  time  for  instabilities  to  grow  within  the 
channel.  Since  most  liquid  rocket  injectors  of  this  type  use  very  small  submergence  lengths,  the 
present  analysis  shows  little  evidence  of  instability  from  the  mechanism  investigated  in  this  work. 

In  general,  the  frequency  of  the  instability  tends  to  decrease  with  increased  submergence.  How¬ 
ever,  the  FFT  of  the  L/D  —  3  case  is  interesting  in  that  two  distinct  frequencies  with  similar 
energy  content  are  noted.  A  short  frequency  consistent  with  a  long-wave  instability  is  present  in 
addition  to  a  higher  frequency  representative  of  clumps  of  disintegrated  fluid.  For  very  short  sub¬ 
mergence  lengths  there  simply  isn’t  time  to  grow  instabilities  to  any  appreciable  level  under  the 
Kelvin-Helmholtz  mechanism. 

Reynolds  number  effects 

The  influence  of  viscous  interactions  were  investigated  through  a  series  of  simulations  at  Reynolds 
numbers  between  103  and  107.  Once  again,  the  other  input  parameters  were  maintained  at  values 
selected  in  the  baseline  configuration.  As  we  can  see  from  the  Figure  17,  the  amplitude  of  the 
oscillation  is  only  slightly  affected  over  this  substantial  Reynolds  number  range.  In  cases  of  high 
Reynolds  numbers,  however,  two  distinct  frequencies  appear  in  the  oscillation;  a  fundamental 
frequency  and  a  first  harmonic  mode  at  a  lower  energy  level.  The  higher  frequency  is  created  by 
liquid  fragments  separated  from  the  liquid  sheet  breakup  near  the  exit.  Since  the  liquid  fragments 
come  out  of  the  nozzle  exit  in  a  very  regular  manner,  the  short  wave  pattern  at  high  frequency  is 
developed.  The  amplitude  of  the  oscillation  at  high  Reynolds  number  decreases  slightly  because 
of  energy  loss  caused  by  the  harmonic  mode  development. 

This  result  is  consistent  with  general  nonlinear  behavior  since  decreased  viscous  interaction 
generates  higher  frequency  mode  in  addition  to  the  primary  oscillation. 
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Effect  of  Liquid  Post  Thickness 

In  this  section  the  “baseline”  flow  conditions  are  considered  with  a  liquid  sheet  with  a  finite  post 
thickness.  This  structure  is  often  referred  to  as  the  “LOX  post”  in  coaxial  injectors  involving 
this  cryogenic  propellant.  When  we  consider  a  finite  thickness  of  the  post,  the  flow  must  now 
expand/recirculate  to  fill  the  region  immediately  downstream  of  the  post  tip.  A  pair  of  counter¬ 
rotating  vortices  become  apparent  immediately  downstream  of  the  upper  and  lower  posts.  The 
recirculation  region  has  approximately  the  same  thickness  as  the  post  tip  thickness  (tiox)  at  the 
nozzle  exit  and  tapers  to  zero  at  some  downstream  point.  This  recirculation  area  serves  as  flame 
holder  in  combusting  system  and  stabilizes  the  liquid  jet  inside  the  injector  by  anchoring  the  jet  to 
the  nozzle  exit. 

Moon  20  investigated  coaxial  injector  and  showed  that  the  recirculation  region  exists  immedi¬ 
ately  after  the  LOX  post  exit.  He  explained  that  the  recirculation  region  is  caused  by  at  least  two 
mechanisms:  viscous  mixing  and  pressure  gradient.  Since  static  pressure  is  not  constant  through¬ 
out  the  flow,  the  recirculation  zone  forms  a  pseudo-diffuser,  causing  a  low  pressure  region  to  exist 
near  the  post  tip.  Viscous  effects  cause  a  slightly  more  rapid  decrease  in  velocity  at  the  jet  bound¬ 
ary. 

Recently,  Glogowski  and  Micci 10  conducted  a  series  of  experiments  and  investigated  the  flow 
near  the  LOX  post  region.  Unfortunately,  measurement  of  the  flow  inside  the  recessed  region  didn’t 
produce  favorable  results  due  to  difficulties  involved  in  characteristically  high  signal-to-noise  ratio 
in  this  area.  However,  they  did  verify  the  recirculation  region  near  the  LOX  post  by  measuring 
mean  axial  velocity  and  showing  that  there  exists  negative  mean  velocity  at  the  post  tip.  For  the 
condition  where  the  LOX  post  thickness  is  not  negligible,  they  agreed  that  macroscopic  instability 
within  the  liquid  jet  arise  from  both  static  pressure  perturbation  and  aerodynamic  viscous  forces. 
In  another  research  of  Glogowski  et  al.  9,  the  flow  characteristics  of  recessed  LOX  post  area  was 
investigated.  They  found  that  the  injection  operation  with  tapered  LOX  post  recessed  into  fuel 
annulus  produced  a  resonance  condition  characterized  by  a  whistling  noise  and  a  significant  mod¬ 
ification  to  the  overall  spray  structure.  The  non-tapered  (straight)  post  exhibited  a  lower  amplitude 
whistling  noise  but  did  not  affect  the  spray  structure. 

All  the  flow  properties  used  in  the  previous  baseline  case  are  preserved  except  to  let  the  liquid 
post  has  a  finite  thickness  with  solid  wall  no-slip  boundary  conditions.  In  order  to  improve  resolu¬ 
tion  in  liquid  post  exit  area,  the  spacing  in  x  direction  is  exponentially  stretched  by  using  the  same 
strategy  used  in  y  direction  stretching.  The  smallest  grid  size  is  applied  to  the  area  between  gas 
and  liquid-phase  where  the  most  interaction  is  expected.  The  liquid  post  thickness  here  is  0.1  each 
in  unit  channel  width.  Dimensional  conversion  gives  each  0.5mm  thickness  of  upper  and  lower 
post  in  5mm  channel  width.  It  should  be  noted  that  no  slip  boundary  condition  is  defined  on  the 
post  tip  in  y  direction.  The  other  boundary  conditions  remain  the  same  as  the  baseline  case  without 
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liquid  post  thickness. 

Figure  18  shows  pseudo-density  contours  at  various  times.  The  liquid  sheet  develops  symmet¬ 
ric  disturbances  on  its  surface  due  to  expansion  of  the  flow  into  the  bluff  region  downstream  of  the 
posts  (a).  The  length  of  this  region  is  about  5  or  6  times  the  post  tip  thickness.  After  this  calcu¬ 
lation  startup  period,  the  liquid  sheet  enters  a  fiber-type  breakup  mode  (Chigier21)  which  creates 
continuous  surface  vortices  showing  aerodynamic  viscous  effects  (b).  This  limit  cycle,  bounded- 
oscillation  state  lasts  until  instability  of  the  liquid  sheet  itself  becomes  prominent  and  starts  to 
generate  asymmetric  wave  (c).  This  wave  grows  rapidly,  and  the  liquid  sheet  oscillates  in  a  quasi- 
periodic  manner  (d),  (e).  When  compared  with  the  baseline  case  without  liquid  post  thickness,  the 
oscillation  is  larger  in  magnitude  and  has  a  shorter  wavelength,  leading  to  a  higher  frequency. 

Velocity  streamlines  and  exit  velocity  profiles  are  shown  in  Figure  19.  The  flow  recirculation 
at  the  liquid  post  exit  is  obvious  in  this  figure.  Symmetric  recirculation  regions  are  developed 
on  upper  and  lower  surfaces  during  the  initial  and  bounded  oscillating  state  period  (a),  (b).  The 
velocity  profiles  at  the  injector  exit  in  this  period  show  the  momentum  deficit  in  the  radial  areas 
corresponding  to  the  upper  and  lower  posts.  During  this  period,  one  can  see  the  assumed  plug-flow 
velocity  distribution  in  the  liquid  is  largely  maintained.  Once  the  liquid  sheet  starts  to  oscillate, 
the  velocity  profile  goes  through  a  corresponding  oscillation.  Note  that  the  highest  density  region 
is  actually  convected  out  of  the  channel  with  the  highest  velocity  in  parts  (d)  and  (e)  even  though 
the  liquid  is  injected  at  lower  velocity  than  the  gas.  These  velocity  profile  characteristics  show  the 
same  qualitative  behavior  as  the  results  of  the  measurements  by  Moon  20  and  Glogowski  et  al.9. 

Two  different  liquid  post  thickness  have  been  evaluated  in  order  to  assess  the  influence  of  this 
parameter.  Posts  of  twice  ( tiox  =  0.2)  and  half  (t lox  =  0.05)  of  the  baseline  thickness  were 
considered.  The  magnitude  of  oscillation  of  each  LOX  post  thickness  are  plotted  with  the  previous 
baseline  case  without  finite  post  thickness  in  Figure  20.  Note  that  the  lines  are  plotted  only  for  the 
oscillation  period,  i.e.,  the  initial  point  (f  =  0.0)  is  when  the  quasi-periodic  oscillation  starts.  For 
the  thicker  post  case,  larger  surface  vortices  are  formed  initially  due  to  bigger  recirculation  region 
at  the  post  tip.  The  liquid  sheet  then  undergoes  a  quasi-periodic  oscillation  with  an  amplitude  much 
less  than  that  of  the  tLox  =  0.1  case.  The  intact  length  of  the  liquid  sheet  from  the  base  of  liquid 
post  is  increased  in  this  case  and  the  frequency  of  the  oscillation  is  lower  than  the  tLox  =  0.1  case. 
In  the  thinner  liquid  post  configuration,  the  long-wave  disturbance  appears  almost  immediately  and 
the  magnitude  of  the  oscillation  grows  to  about  the  same  level  as  the  tLOX  =  0.1  case,  but  with  a 
significantly  higher  frequency  of  oscillation. 

In  comparing  these  results  to  the  case  from  the  prior  section  where  the  post  thickness  was 
neglected  altogether,  we  see  an  interesting  behavior.  The  thicker  post  shows  a  stabilizing  effect, 
while  the  two  thinner  posts  show  amplification  in  the  instability  as  compared  to  the  case  where 
a  zero  post  thickness  was  assumed.  These  results  indicate  that  that  the  thickness  of  the  post  can 
play  a  strong  role  in  the  instability  and  can  in  principle  explain  the  significant  differences  in  spray 


111 


character  noted  by  Glogowski,  et.  al.  when  tapered  and  untapered  posts  were  used. 


Conclusions 

A  homogeneous  flow  model  has  been  developed  to  assess  hydrodynamic  instabilities  of  coaxial  at¬ 
omizers  in  which  a  liquid  jet/sheet  is  submerged  slightly  in  an  annular  gas  stream.  The  calculations 
have  been  performed  in  a  2-D  analogue  to  this  geometry;  a  liquid  sheet  with  coflowing  gas  on  upper 
and  lower  extremities.  Kelvin-Helmholtz  type  instabilities  are  noted  under  conditions  where  the 
gas  velocity  is  substantially  greater  than  that  of  the  liquid.  Parametric  studies  indicate  the  ampli¬ 
tude  of  the  instability  increases  with  gas  velocity,  gas  density,  and  sheet  recess/submergence  length 
inside  the  channel.  Increasing  sheet  thickness  tended  to  decrease  the  amplitude  of  the  oscillation. 
The  frequencies  observed  are  consistent  with  the  channel  transit  time  of  the  liquid. 

A  series  of  liquid  post  thicknesses  were  investigated  to  assess  the  effects  of  this  parameter. 
Flow  recirculation  at  the  post  tip  is  observed  in  accordance  with  prior  experimental  observations. 
Thin  posts  show  greater  amplitude  instabilities  than  the  zero-thickness  case,  while  a  thicker  post 
proved  to  damp  the  amplitude  of  the  oscillation.  In  addition,  the  thinner  posts  created  much  higher 
oscillation  frequencies  than  the  case  of  a  zero-thickness  post.  These  results  are  indicative  of  the  ef¬ 
fects  noted  experimentally;  strong  differences  in  exit  flow  can  be  attributed  to  fairly  minor  changes 
in  post  thickness. 
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Nomenclature 

D  -  Channel  width 

/  -  Nondimensional  frequency 

h  -  Liquid  sheet  thickness  from  center  line 

k  -  Wave  number  L  -  Channel  length 

m  -  Nozzle  mass  flow  rate 

P  -  Pressure 

Re  -  Reynolds  number 

tLOx  -  Liquid  post  thickness 

u,  v  -  Velocity 
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V  -  Computational  cell  volume 
We-  Weber  number 
x  -  Cartesian  coordinate,  x-direction 
y  -  Cartesian  coordinate,  y-direction 
a  -  Void  fraction 

A  U  -  Velocity  difference  between  liquid  and  gas  phase 
c  -  Density  ratio  ( pg/pi ) 

77  -  amplitude  of  wave  disturbance 

H  -  Viscosity 

<t>  -  Velocity  potential 

p  -  Fluid  pseudo-density 

a  -  Surface  tension 

u  -  Growth  rate 


Subscripts 

c  -  Computational  cell 
g  -  Gas-phase 
l  -  Liquid-phase 
i  -  Inlet 
o  -  Outlet 
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Figure  Captions 

1.  Nomenclature  for  linear  stability  analysis  for  anti-symmetric  jet  oscillations  in  a  channel 


2.  Effect  of  gas  velocity  on  growth  rate:  e  =  0.01,n  =  5 

3.  Effect  of  gas/liquid  density  ratio,  e  on  growth  rate:  A U  =  3,  n  =  5 

4.  Effect  of  channel  width  on  growth  rate:  e  =  0.01,  A U  =  3 

5.  Coaxial  injector  schematic  noting  design  variables. 

6.  Effect  of  grid  size  on  density,  pressure,  and  velocity  at  a  point  on  the  centerline  of  the  domain 
(Ug/Ui  =  4,  L/D  =  2,  h  =  0.2,  Ret  =  105) 

7.  Density  distribution  at  exit  at  steady  state 

8.  Velocity  distribution  at  exit  at  steady  state 

9.  Density  contours  in  a  typical  jet  oscillation 

10.  Velocity  streamlines  and  exit  velocity  profiles  in  a  typical  jet  oscillation 

11.  Liquid  sheet  oscillation  :  Maximum  density  location  at  channel  exit  ( Ug/Ui  =  4,  L/D  = 
2,  h  =  0.2,  Rei  =  105) 

12.  Fast  Fourier  Transform  analysis  of  liquid  sheet  oscillation 

13.  Effect  of  density  ratio  on  amplitude  of  the  sheet  oscillation  ( Ug/Ui  =  4,  L/D  =  2,  h  = 
0.2,  Ret  =  105) 

14.  Effect  of  velocity  ratio  on  amplitude  of  the  sheet  oscillation  (e  =  0.01,  L/D  =  2,  h  = 

0.2,  Rei  =  105) 

15.  Effect  of  jet  thickness  on  amplitude  of  the  sheet  oscillation  ( Ug/Ui  =  4,  e  =  0.01,  L/D  = 

2,  Ret  =  105) 

16.  Effect  of  channel  length  on  amplitude  of  the  sheet  oscillation  ( Ug/Ui  =  4,  e  =  0.01,  h  = 

0.2,  Rei  =  105) 

17.  Effect  of  gas  Reynolds  number  on  amplitude  of  the  sheet  oscillation  ( Ug/Ui  =  4,  e  =  0.01, 
L/D=2,  h  =  0.2) 

18.  Density  contours  liquid  sheet  oscillation  with  liquid  post  thickness,  tLox  =  0.1 
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19.  Velocity  streamlines  and  exit  velocity  profiles  with  liquid  post  thickness,  tLOX  =  0.1 

20.  Effect  of  liquid  post  thickness  on  amplitude  of  oscillation  (Ug/Ui  =  4,  L/D  =  2) 
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Figure  3:  Effect  of  gas/liquid  density  ratio,  e  on  growth  rate:  AU  =  3,  n 


Figure  4:  Effect  of  channel  width  on  growth  rate:  e  =  0.01,  AU  =  3 
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Computational  Domain 


Figure  5:  Coaxial  injector  schematic  noting  design  variables 


a)  Density  history 


b)  Pressure  history 


c)  Velocity  history 


Figure  6:  Effect  of  grid  size  on  density,  pressure,  and  velocity  at  a  point  on  the  centerline  of  the 
domain  ( Ug/Ui  =  4,  L/D  —  2,  h  =  0.2,  Ret  =  105) 


Figure  7:  Density  distribution  at  exit  at  steady  state 


Figure  8:  Velocity  distribution  at  exit  at  steady  state 
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Figure  10:  Velocity  streamlines  and  exit  velocity  profiles  in  a  typical  jet  oscillation 
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Figure  11:  Liquid  sheet  oscillation  :  Maximum  density  location  at  channel  exit  ( Ug/Ui 
4,  L/D  =  2,  h  =  0.2,  Ret  =  105) 


Figure  12:  Fast  Fourier  Transform  analysis  of  liquid  sheet  oscillation 
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Figure  13:  Effect  of  density  ratio  on  amplitude  of  the  sheet  oscillation  (Ug/Ui  =  4,  L/D  =  2,  h 
0.2,  Ret  =  105) 


Figure  14:  Effect  of  velocity  ratio  on  amplitude  of  the  sheet  oscillation  (e  =  0.01,  L/D  =  2,  h 
0.2,  Ret  =  105) 
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Figure  15:  Effect  of  jet  thickness  on  amplitude  of  the  sheet  oscillation  (Ug/Ui  =  4,  e  = 
0.01,  L/D  =  2,  Ret  =  105) 


Figure  16:  Effect  of  channel  length  on  amplitude  of  the  sheet  oscillation  ( Ug/Ui  =  4,  e  = 
0.01,  h  =  0.2,  Ret  =  105) 


Figure  18:  Density  contours  liquid  sheet  oscillation  with  liquid  post  thickness,  t^ox  —  0.1 
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9  Appendix  E  -  Three-Dimensional  Coaxial  Injector  Modeling 


Kim*,  B.  and  Heister,  S.  D.,  “Three-dimensional  Simulations  of  Flow  within  the  Recessed 
Region  in  a  Coaxial  Injector”. 
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Abstract 

A  three-dimensional,  two-phase,  unsteady  Navier-Stokes  solver  has  been  developed  to  investigate 
fluid  dynamic  instabilities  within  the  recessed  region  of  a  shear  coaxial  injector  element.  Here,  the 
main  emphasis  is  to  study  applications  related  to  liquid  rocket  engine  injectors  using  the  gas/liquid 
shear-coaxial  element  in  which  the  inner  liquid  cylindrical  post  is  submerged  slightly  with  respect 
to  the  overall  exit  plane  of  the  device.  Since  most  of  the  previous  works  focus  on  spray  structure 
outside  of  the  injector,  this  study  provides  readers  with  insight  into  unsteadiness  resulting  from 
hydrodynamic  instabilities  within  the  internal  nozzle  flow  upstream  of  the  combustion  chamber. 

The  present  study  focuses  on  unsteady  ’self-oscillations’  which  have  been  theorized  by  various 
researchers.  The  Kelvin-Helmholtz  instability  mechanism  due  to  velocity  discontinuity  between 
gas  and  liquid  phase  is  investigated  as  a  source  of  unsteadiness  which  could  contribute  to  combus¬ 
tion  roughness  or  instabilities.  Massflow  variations  of  the  order  of  30-40%  are  shown  to  exist  as  a 
result  of  this  highly  nonlinear  process;  fundamental  frequencies  are  also  identified  for  a  range  of 
conditions. 


Introduction 

The  liquid-oxygen/gaseous-hydrogen  coaxial  injectors  have  been  widely  used  in  the  liquid  rocket 
engines.  The  concentric-orifice  coaxial  injector  was  used  in  rocket  engine  testing  as  early  as  the 
1940s  and  later  became  the  general  choice  on  most  of  cryogenic  liquid  rocket  engines.  The  mech¬ 
anism  of  typical  shear  coaxial  injector  is  that  the  liquid  oxydizer  is  fed  through  the  central  tube 
and  the  gaseous  hydrogen  is  coflowing  around  the  outer  annulus.  Often  times,  the  liquid  oxygen 
post  structure  is  recessed  from  the  injector  face  in  order  to  improve  the  combustion  stability.  Inner 
liquid  turbulence  and  gas  to  liquid  interactions  cause  the  liquid  to  be  stripped  away  from  the  jet 
and  entrained  into  the  surrounding  gasflow.  Disturbances  from  inside  the  combustion  chamber 
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or  outside  vibrations  during  engine  operations,  along  wioth  the  intrinsic  jet  instability  inside  the 
nozzle,  create  a  certain  circumstance  that  the  liquid  jet  cannot  maintain  a  uniform  distribution  of 
atomized  spray.  The  spray  jet  in  the  recessed  region  then  could  produce  periodic  fluctuations  in  the 
combustion  zone  which  may  or  may  not  couple  to  the  dynamics  within  the  combustion  chamber, 
possibly  followed  by  severe  combustion  instability. 

In  general,  for  the  flow  at  high  velocity,  most  of  researchers  agree  that  the  principal  source 
introducing  instability  to  the  jet  is  from  aerodynamic  forces  arising  from  the  interaction  of  the 
liquid  jet  with  the  surrounding  gas  flow.  Reynolds  and  Weber  numbers  are  generally  very  high 
in  these  atomizers  and  aerodynamic  forces  are  several  orders  of  magnitude  larger  than  capillary 
forces.  The  interaction  between  the  liquid  and  gas  phases  mainly  comes  from  different  velocities 
of  each  phase.  The  velocity  discontinuity  in  a  homogeneous  fluid  results  wave  growth  on  the 
interface,  which  is  a  common  Kelvin-Helmholtz  instability. 

The  historical  review  of  the  development  of  various  coaxial  injectors  used  in  liquid  engines 
can  be  found  in  Ref.  1  and  2.  Theoretical  and  numerical  modeling  efforts  on  the  coaxial  injector 
atomization  is  also  summarized  in  the  paper  of  Vingert  et  al.3  along  with  experimental  results. 
Their  work  mainly  devoted  to  the  drop  size  prediction  by  using  those  experimental  data.  In  the 
paper  of  Hulka  and  Hutt 1  the  review  focused  on  instability  phenomena  associated  with  combustion 
chamber  feed  system.  Series  of  injector  testing  revealed  that  injection  velocity  ratio,  pressure 
drop  after  the  injection  and  fuel  injection  temperature  are  the  key  parameters  defining  combustion 
instability. 

According  to  the  Ref.  2,  the  combustion  instability  with  high  frequency  is  categorized  into 
acoustic  instability  and  hybrid  instability.  The  acoustic  instability  shows  dominant  wave-type 
oscillation  in  the  main  chamber,  but  is  independent  of  the  feed  system.  With  the  hybrid  form 
of  instability,  the  wave  character  of  the  oscillation  is  strongly  coupled  between  the  feed  system 
and  the  combustion  chamber.  Hutt  and  Rocker  4  also  investigated  the  high  frequency  combus¬ 
tion  instability  associated  with  coaxial  injectors.  They  classified  the  instability  phenomena  in  the 
chamber  as  injection-coupled  and  intrinsic  mechanism.  The  injection  coupling  implies  chamber 
pressure/temperature  variation  as  a  key  contributor  in  the  change  of  flow  dynamics  through  the 
injector.  In  the  other  hand,  the  intrinsic  mechanism  occurs  in  the  flowfield  due  to  its  own  flow 
dynamics  with  negligible  feed  system  effect.  However,  a  general  agreement  does  not  exist  in  term 
of  categorizing  those  instability  mechanism  due  to  the  involvement  of  various  subprocesses  occur¬ 
ring  after  the  injection  at  different  time  scales.  It  should  be  noted  that  injection  coupling  is  never 
independent  of  the  intrinsic  subprocesses,  such  as  atomization,  propellant  heatup,  vaporization, 
and  mixing  because  these  processes  determine  the  relationship  between  the  injector  response  and 
the  chamber  response. 

Mayer  and  his  research  group  have  done  a  significant  amount  of  work  on  the  coaxial  injector 
in  terms  of  the  combustion  instability.  Mayer  and  Krulle  5  investigated  coaxial  flow  mixing  phe- 
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nomena  in  terms  of  chamber  pressure  variation,  density/velocity  ratio  changes  and  surface  tension 
effect.  By  increasing  chamber  pressure  gas  density  is  increased,  magnifying  the  aerodynamic  in¬ 
teraction  between  the  liquid  and  gas  phases  and  resulting  in  a  faster,  finer  atomization.  Increasing 
surrounding  gas  velocity  also  leads  to  an  increase  of  surface  wave  growth  and  to  macroscopic  in¬ 
stabilities  of  the  liquid  jet.  They  claimed  the  initiation  of  the  jet  surface  deformation  was  due  to 
internal  liquid  turbulence  delivering  energy  transform  in  forms  of  eddy  structures,  approximately 
a  size  of  10-30%  of  the  LOX  post  diameter.  Their  global  spray  simulation  combining  primary  and 
secondary  jet  breakup  showed  qualitative  match  with  their  experimental  results.  In  other  works 
of  Mayer  6,r,  the  coaxial  injector  flow  was  studied  under  firing  engine  condition  at  supercritical 
chamber  pressure  (higher  than  5  MPa).  The  study  revealed  a  remarkable  difference  between  sub- 
critical  spray  formation  and  the  supercritical  injection  and  mixing.  At  subcritical  condition,  the 
liquid  jet  shows  similar  flow  pattern  to  the  cold  flow  test  forming  ligaments  off  the  liquid  jet  sur¬ 
face  and  producing  droplets  before  evaporation.  Upon  approaching  and  exceeding  supercritical 
pressure,  droplets  no  more  exist  and  the  liquid  jet  rapidly  dissolves.  The  flame  from  combustion 
chamber  was  anchored  at  the  tip  of  LOX  post  by  flow  recirculation  eddies  serving  as  flame  holder 
for  steady-state  combustion.  The  asymmetric  flow  oscillation  was  also  reported  in  all  experiments, 
but  the  source  of  the  oscillation  was  not  clearly  stated. 

For  analysis  of  the  wave-type,  asymmetric  jet  oscillation  many  researchers  studied  its  theo¬ 
retical  aspects  and  conducted  numerous  experiments.  Eroglu  and  Chigier  9  focused  on  the  wave 
characteristics  of  liquid  jet  from  coaxial  air-blast  injector.  They  measured  frequency  and  wave¬ 
length  of  the  jet  issuing  from  the  injector,  and  found  two  dominant  wave  type;  spanwise(dilational) 
and  stream-wise(sinuous)  waves.  The  span  wise  wave  showed  at  a  low  relative  jet  velocity  between 
gas  and  liquid-phase  and  the  stream-wise  wave  showed  at  a  high  relative  jet  velocity.  Average 
wavelengths  decreased  with  liquid  and  gas  velocity.  The  frequency  band  of  the  jet  oscillation  in¬ 
creased  with  the  liquid  jet  velocity.  Mansour  and  Chigier 10  also  conducted  similar  research  on  the 
liquid  sheet  instability  issuing  from  the  two-dimensional  air-assisted  nozzle.  The  results  showed 
the  same  pattern  with  the  Eroglu  and  Chigier’s  study.  However,  Mansour  and  Chigier’s  study 
probed  higher  velocity  cases,  and  confirmed  that  the  frequency  of  the  liquid  sheet  oscillation  in¬ 
creased  with  coflowing  gas  velocity.  This  also  indicates  that  the  aerodynamic  interaction  between 
gas  and  liquid  flow  is  the  dominant  factor  for  the  flow  oscillation. 

Instability  mechanisms  in  coaxial  injectors  were  investigated  experimentally  by  Glogowski  et 
al.  13-15  under  noncombusting  conditions.  Their  experiment  results  showed  that  for  the  coaxial 
injector  with  the  liquid  oxygen  (LOX)  post  recessed  into  the  fuel  annulus,  the  injector  transitioned 
into  a  condition  of  resonance  characterized  by  a  whistling  noise  and  significant  modification  to  the 
overall  structure  of  the  spray  due  to  the  strong  acoustic  coupling  between  injector  hydrodynamics 
and  spray  formation.  Without  the  recessed  region,  the  injector  operation  produced  a  near  resonance 
condition  with  a  lower  amplitude  whistling  noise  but  did  not  make  considerable  changes  in  the 
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spray  structure.  The  recirculation  by  reverse  flow  near  the  recessed  LOX  post  exit  also  confirmed 
in  the  study.  A  try  for  getting  oscillation  frequency  data  in  the  submergence  region  of  the  coaxial 
injector  was  not  successful  due  to  optical  restraint  with  small  scales  involved. 

Bazarov  16-18  investigated  coaxial  injector  flow  dynamics  in  engine  operating  conditions.  He 
identified  “self-oscillation”  and  “self-pulsation”  modes  of  liquid  jet  characteristics  and  suggested 
them  as  sources  of  high  amplitude  noise  during  combustion,  leading  to  combustion  inefficiency. 
He  claimed  that  those  high  frequency  instability  occurred  because  of  the  intrinsic  unsteadiness 
of  flowfield  and  interactions  with  the  combustion  chamber  and  feed  system  dynamics  since  the 
injector  element  operates  in  strong  feedback  loops  coupling  chamber  response  and  injection  feed 
system.  However,  it  was  not  clear  which  mechanism  had  more  serious  influence  on  the  jet  oscilla¬ 
tion/pulsation. 

Most  of  the  previous  research  works  summarized  here  mainly  contributed  to  the  investigation 
of  flow  phenomena  outside  of  injector  after  the  injection.  None  of  them  addressed  the  inner  flow 
structure  in  a  recessed  region  of  the  coaxial  injector,  especially  using  numerical  simulation.  Since 
the  upstream  flow  in  a  injector  provides  the  initial  condition  for  the  entire  spray  atomization  process 
and  combustion  chamber  acoustic  characteristics,  research  work  on  this  area  is  highly  desirable. 

Motivated  by  Bazarov’s  research,  Kim  and  Heister  19-21  developed  a  2-D  numerical  model 
investigating  the  coaxial  liquid  sheet  behavior.  The  simulation  demonstrated  the  self-oscillation 
mode  caused  by  the  Kelvin-Helmholtz  instability.  A  series  of  parametric  study  had  been  done 
in  order  to  analyze  the  effects  of  several  key  parameters  including  velocity/density  ratio,  recess 
length  changes  and  the  LOX  post  thickness.  The  results  showed  that  increasing  density  and  ve¬ 
locity  ratio  produced  higher  oscillation  frequency  and  increasing  the  recess  length  also  increased 
the  natural  frequency  of  the  liquid  sheet  oscillation,  leading  to  the  sheet  breakup.  These  results 
showed  the  same  tendency  that  was  predicted  by  linear  theory  analysis.  Recent  experimental  study 
of  Branam  22  provided  the  same  qualitative  match  to  the  overall  liquid  jet  behavior.  Based  on  the 
2-D  results,  a  three-dmensional  model  has  also  developed  and  a  series  of  has  been  conducted  using 
Space  Shuttle  Main  Engine’s  injector  geometry  and  its  flow  properties.  In  this  paper,  the  hydrody¬ 
namic  instability  of  the  coaxial  jet  in  the  recessed  region  is  presented  using  the  three  dimensional 
direct  numerical  simulation.  The  results  indicate  that  the  recess  length  to  injector  orifice  diameter 
has  a  significant  effect  on  spray  structure  over  the  velocity  and  gas  ratio  changes.  Unfortunately, 
none  of  experimental  data  for  exact  comparison  exists  due  to  very  small  spartial  and  temporal 
scales  involved  in  the  area.  However,  some  of  recent  experimental  data  provide  good  measure  in  a 
macroscopic  point  of  view  for  the  instability  frequency  as  explained  later  in  this  paper. 
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Modeling  Description 


Numerics 

A  fully  unsteady,  three  dimensional,  two-phase  simulation  has  been  developed  utilizing  a  finite 
volume  implementation  of  the  Marker  and  Cell  discretization  method.  The  current  model  is  based 
on  locally  homogeneous  flow  (LHF)  assumption  in  which  the  relative  velocity  and  temperature 
between  two-phases  are  small  enough  in  comparison  to  variation  of  the  overall  flow  field  that  is 
predicted.  Additional  constituitive  relation  for  desity  field  has  been  implemented  in  order  to  pro¬ 
vide  a  mechanism  of  solving  two-phase  flow  with  a  single  phase  Navier-Stokes  equations  set.  This 
fictious  “pseudo”  density  varies  in  amplitude  between  the  liquid  and  gas  extream.  The  LHF  as¬ 
sumption  and  the  pseudo-density  implementation  allow  the  current  model  to  handle  the  two-phase 
flow  field  with  one  governing  equations  set  rather  than  to  compute  separate  governing  equation 
sets  for  each  flow  phase,  liquid  and  gas  in  this  case. 

Formulation  of  the  three  dimensional  governing  equations  for  computation  are  given  as  follows. 


dpw  dpuw  dpvw  dpw 2  dP 
dt  +  dx  +  dy  +  dz  +  dy 


d  dw  d  dw  d  dw 
dx^  dx  dy^  dy  dz^  dz 


The  Lagrangian  form  of  the  continuity  equation  is  also  required  : 

Dp  ,du  dv  dw .  ... 

TTt+*Tx  +  irv  +  aI)  =  0  (5) 

Under  the  two-phase  flow  condition,  the  viscosity  can  vary  spatially.  According  to  Kubota  et  al. 
26 ,  the  viscosity  of  mixture  can  be  written: 

p  =  apg  +  (1  -  a)pi  (6) 

where  pg  and  pi  are  the  gas  and  liquid  viscosities,  and  a  is  the  void  fraction.  Since  the  non- 
dimensional  pseudo-density  is  volume  fraction  of  mass  per  unit  cell  volume,  the  Equation  (6)  can 
be  written  as: 

p(p)  =  pm  +  (1  -  p)Pg  (7) 
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This  mixture  viscosity  is  substituted  back  into  Equation  (2)  through  (4)  for  non-dimensionalization. 
The  channel  width,  liquid  inflow  velocity  and  liquid  density  are  chosen  as  dimensions  in  nondi- 
mensionalizing  the  equations.  Then  the  momentum  equations  (2  -  4)  can  be  arranged  as  follows 
for  the  momentum  flux  calculation. 
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where  the  Reynolds  number  for  mixture  Re*  is  defined  as  follow: 
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which  includes  the  both  liquid  and  gas  phase  viscous  effect  in  one  term.  Here  Rei  and  Reg  repre¬ 
sent  Reynolds  numbers  of  liquid  and  gas-phase  respectively. 

An  additional  constituitive  relation  for  the  non-physical  pseudo-density  variable  is  required  in 
place  of  an  equation  of  state  in  order  to  close  the  set  of  governing  equations.  Basically,  the  droplets 
in  gas  media  are  assumed  to  be  the  same  size,  i.e.  the  droplet  volume  does  not  change  due  to  the 
pressure  variation  of  the  flow  field.  Coalescence  or  collision  between  the  droplets  are  neglected, 
nor  the  secondary  breakup  of  droplet  is  considered.  Consequently,  one  can  conclude  from  the 
above  assumptions  that  the  density  evolves  in  time  by  simply  considering  a  Lagrangian  tracking  of 
the  droplet  field  as  specified  by  a  continuity  equation  as: 


Dp 

Dt 


dp 

dt 


dp  dp  dp 
+  u-£-  +  v-£-  +  w—  =  0 
dx  dy  dz 


(12) 


This  equation  is  basically  a  statement  that  the  droplet  mass  is  invariant  along  a  path  line  in  the 
fluid.  By  taking  account  for  the  flow  direction,  the  pseudo-density  is  updated  based  on  mass  flux 
calculation  as  follows : 
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where  Vc  is  volume  of  a  given  cell  and  rh  is  mass  flowrate  corresponding  x,  y  and  z  directions. 
The  flux  of  the  mass  through  a  given  cell  are  calculated  by  a  standard  upwind  scheme.  While 
this  provides  a  locally  first  order  solution  in  the  region  adjacent  to  the  discontinuity  in  density,  it 
provides  for  a  stable  approach  to  account  for  the  large  density  gradients  near  the  interface.  Since  the 
interface  region  will  undoubtedly  contain  droplets  in  these  very  high  convective  environments,  a 
shock-type  density  discontinuity  would  not  be  consistent  with  a  homogeneous  flow  representation. 
The  scheme  is  second-order  accurate  with  the  exception  of  the  points  adjacent  to  the  interface, 
and  convergence  is  verified  in  the  next  section.  One  can  find  more  details  on  development  of  the 
numerical  solver  in  Bunnell’s  work  25 .  In  this  study,  the  3-D  coaxial  jet  modeling  work  will  be 
presented  focusing  on  hydrodynamic  instability  of  shear  coaxial  jet. 

Parallelization 

The  code  runs  on  a  state-of-art  Beowulf  Linux  cluster  that  is  equipped  with  104  processors  and  fast 
ethemet  network.  One  run  usually  requires  12  to  24  processors  depending  on  mesh  discretization. 
Since  the  cluster  is  dedicated  to  the  modeling,  each  run  makes  use  of  nearly  100%  of  CPU  power 
and  network  bandwidth.  Even  with  this  superb  environment,  one  run  up  to  150,000  time  steps 
takes  about  three  weeks. 

Parallel  processing  using  MPI  (Message  Passing  Interface)  has  been  implemented  in  order  to 
run  the  3-D  model  in  a  timely  manner.  Even  with  the  large  number  of  1 .2GHz  AMD  Athlon 
processors  used  in  the  study,  running  a  full  3-D  CFD  model  in  a  timely,  economic  manner  still 
remains  as  a  challenge.  The  MPI  has  become  a  standard  method  for  the  parallel  processing  along 
with  PVM  (Parallel  Virtual  Machine).  For  this  study,  the  MPI  has  been  adopted  due  to  its  easy-to- 
use  and  built-in  libraries  in  many  standard  Linux  distributions  such  as  RedHat. 

MPI  libraries  provide  two  basic  routines,  a  function  for  sending  data  and  one  for  receiving, 
along  with  a  multitude  of  other  useful  functions  that  provide  the  collective  communication  opera¬ 
tions  that  are  necessary  for  a  group  of  processors  running  parallel  processing  jobs. 

For  the  present  3-D  atomization  modeling,  the  computational  domain  is  split  up  in  axial  direc¬ 
tion  for  the  desired  number  of  processors,  n.  While  each  processor  solves  a  flow  field  of  subdivided 
domain,  the  boundary  conditions  of  each  subdomain  are  transfered  to  neighboring  domain  through 
message  passing.  All  the  processors  involved  in  this  procedure  perform  the  same  calculation  in 
each  step  by  copying  the  original  task  to  other  processors  but  solving  different  subdomains  and 
boundary  conditions.  This  type  of  parallelism  is  called  SPMD,  which  stands  for  Single-Program- 
Multiple-Data.  A  typical  example  of  using  the  SPMD  parallel  process  is  to  solve  the  Poisson 
equation. 

In  the  finite  volume  Marker  and  Cell  method  used  in  this  study,  the  mesh  discretization  is  such 
that  the  calculation  of  primitive  variables,  P,  p,  V,  and  the  Lagrangian  derivative  Dp/Dt,  only 
depend  on  values  at  neighboring  nodes.  Thus,  it  is  required  that  boundary  values  on  a  neighboring 
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node  be  known  for  the  calculation  of  primitive  variables,  and  the  data  along  the  boundary  of  a 
given  domain  must  be  passed  to  its  neighbor.  Physically,  the  problem  being  solved  in  parallel  on 
V  processors  can  be  viewed  as  V  separate  flow  field  problems  each  being  solved  with  different 
set  of  boundary  conditions.  However,  the  boundary  conditions  for  each  problem  are  dependent  on 
the  flow  field  in  each  domain.  For  processor  CP\  the  inlet  boundary  condition  is  a  constant  flow 
field  condition,  as  is  the  exit  boundary  condition  for  the  flow  field  being  solved  by  CPn_i.  Here, 
the  CPi  in  this  section  denotes  the  computing  processor  corresponding  to  each  subdomain  to  be 
solved.  The  inner  domain  boundary  conditions  are  function  of  time  at  the  inlet  and  exit  and  are 
functions  of  the  surrounding  domain  flow  field. 

Figure  1  shows  the  schematic  of  computational  domain  discretization.  Assuming  that  com¬ 
putational  domain  has  m2  x  n2  grid  points  in  i-  and  j-  directions,  which  represent  x  and  y 
coordinates  respectively  in  physical  domain  along  with  k-  for  z  direction,  and  that  “m2”  in  ax¬ 
ial  ( i -)  direction  is  evenly  divisible  by  number  of  processor  “n”,then  each  domain  will  contain 
(m2/n)*n2  grid  points  in  the  streamwise  direction.  Since  the  boundary  values  at  the  sides  of  a 
given  domain  are  needed  by  other  processors  due  to  differencing  scheme  used,  the  data  size  for 
one  primitive  value  in  each  domain  will  be  ((m2/n)  +  2)  *  n2.  Here,  the  initial  and  last  data  array 
where  i  =  1  and  i  =  (m2/n)  +  2  would  be  locations  for  the  boundary  data  that  is  passed  between 
processors,  and  the  arrays  between  those  two  locations  (from  i  =  2  to  i  =  (m2/n)  +  1)  would  be 
the  domain  to  be  solved.  A  3-D  problem  requires  the  passing  of  arrays  for  the  boundary  condition 
transmission  while  2-D  problem  requires  the  passing  of  vectors.  The  total  data  size  required  for 
all  processes  would  be  n  *  ((m2/n)  +  2)  *n2*  nbyte  bytes,  which  depending  on  whether  the 
variables  are  declared  as  single  precision  or  double  precision,  nbytes  would  be  equal  to  2  or  4 
bytes  respectively. 

Another  issue  associated  with  parallelization  in  this  modeling  is  the  Poisson  solver  for  calcu¬ 
lation  of  pressure  field.  The  system  of  sparse  linear  equations  has  to  be  solved  in  such  a  way  that 
the  solution  satisfies  each  process  domain.  The  original  solution  algorithm  in  2-D  simulations 19,20 
relies  on  a  Successive  Line  Over  Relaxation  (SLOR)  method,  and  the  parallelized  solver  takes  the 
same  methodology.  The  SLOR  technique  is  applied  to  each  subdomain  independently,  and  yet 
each  subdomain  is  linked  to  others  through  the  message  passing.  The  procedure  is  outlined  as 
follows  referring  to  Figure  1; 

•  Sweep  the  first  column  in  each  domain. 

•  Pass  the  result  from  CPn+l  to  CPn. 

•  Finish  the  sweep  across  the  domain. 

•  Pass  the  result  from  CPn  to  CPn+\. 
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This  provides  linkage  of  the  domains  at  the  boundaries  since  the  SLOR  sweep  requires  knowing 
the  value  of  P  at  locations  “z  +  1”  and  “z  —  1”  of  each  domain.  The  communication  requirements 
are  such  that  two  calls  to  MPI  subroutines,  MPI.SEND  and  MPIJRECV,  must  be  made  for  each 
iteration.  References  28  -  32  provide  background  information  on  MPI. 

Despite  the  large  number  of  processors  available,  speed  was  maximized  when  6  grid  points 
(minimum  number  of  grid  points  decided  due  to  differential  schemes  used  in  the  solver)  in  axial 
direction  for  each  individual  processors  are  used,  resulting  12  processors  of  use  for  72  x  71  x  71 
mesh.  With  the  individual  processors  at  this  speed  (1.2  GHz  chips),  the  computation  becomes 
communcation  limited  when  larger  number  of  nodes  are  employed. 

Grid  Refinement  Tests 

The  computational  domain  consists  of  the  liquid  oxydizer  (LOX)  and  gaseous  hydrogen  flow  path 
from  the  LOX  post  exit  to  the  injector  orifice  exit,  which  is  the  combustion  chamber  face.  To 
account  for  the  variation  in  the  flowfield,  the  domainu  is  devided  into  numerous  grid  cells  such  that 
the  volume  of  a  computational  cell  is  much  smaller  than  the  wavelength  of  the  highest  frequency 
acoustic  mode  to  be  considered. 

The  computational  structured  mesh  generation  consists  of  two  parts;  one  is  inner  circle  for 
liquid-phase,  and  the  other  is  outer  annulus  for  gas-phase  flow.  Due  to  the  liquid  post  structure 
configured  between  the  liquid  cylinder  and  gaseous  fuel  annular  passage,  the  circular  shape  of 
physical  domain  has  to  be  fixed.  However,  a  polar  coordinate  mesh  cannot  be  used  due  to  the 
condition  of  maintaining  n  x  n  sparse  matrix  solver.  The  inner  circle  meshing  begins  with  the 
meshing  of  a  square.  The  square  mesh  is  then  deformed  making  use  of  polynomial  stretching  of 
order  n  in  order  to  form  a  quarter  of  a  circle.  A  mirror  image  is  then  symmetrically  projected  until  it 
forms  a  complete  circle.  The  reason  that  we  do  not  extend  this  out  to  the  nozzle  wall  is  because  we 
have  to  have  complete  the  circle  configuration  for  the  liquid  post  structure.  For  the  outer  gas-phase 
domain,  circles  surround  the  inner  liquid  part  while  lines  in  y  and  z  directions  are  extended  from 
the  nodes  on  outest  circle  of  inner  domain  to  the  nozzle  wall  side.  A  complete  cross-sectional  view 
of  a  coarse  computational  mesh  is  shown  in  Figure  2.  Grid  points  in  the  diagonal  directions  create 
challenges  since  the  n  x  n  matrix  size  should  be  kept  due  to  the  flow  solver  algorithm.  This  problem 
was  solved  by  deploying  two  lines  in  both  y  and  2  directions  from  the  diagonal  nodes  to  external 
circumference,  resulting  nearly  triangular-shaped  cells  on  the  diagonal  directions.  The  liquid  post 
tip  area  has  the  smallest  grid  size  to  resolve  the  recirculation  zone  immediately  downstream  of  the 
post. 

A  grid  convergence  test  has  been  conducted  in  order  to  find  an  adequate  mesh  size.  Three  mesh 
sizes:  30  x  30  x  30,  60  x  60  x  60  and  90  x  90  x  90  in  x,  y  and  z  direction  respectively  are  tested 
and  the  result  is  shown  in  Figure  3.  However,  due  to  the  run-time  constraint,  the  test  is  limited 
to  a  short  period  of  time  up  to  10,000  time  steps.  Figure  3  depicts  the  time  history  of  mass  flow 
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through  the  injector  nozzle  exit  for  each  mesh  size.  While  the  coarse  grid  (30  x  30  x  30)  shows 
unacceptable  behavior  with  too  much  of  diffusion,  the  middle  size  grid  (60  x  60  x  60)  matches 
fairly  well  with  the  fine  grid  (90  x  90  x  90).  Though  running  with  the  fine  grid  or  even  finer  grid 
would  be  desirable,  the  cost  for  running  those  fine  grids  requires  unbearable  run-time  (more  than 
one  month  up  to  10,000  time  steps  for  90  x  90  x  90  grid  size).  Therefore,  the  middle  grid  size  has 
been  adopted  for  the  3-D  modeling.  Actual  3-D  simulations  have  used  a  72  x  71  x  71  mesh  by 
increasing  the  number  of  grid  points  as  many  as  possible  without  sacrificing  a  reasonable  run-time. 

Boundary  conditions  are  applied  to  the  domain  as  follows:  Along  the  wall  of  the  injector, 
the  no-slip  condition  is  used.  The  liquid  post  tip  area  at  the  inlet  plane  also  requires  a  no-slip 
condition.  On  the  inflow  boundary,  constant  liquid/gas-phase  velocities  are  defined.  The  pressure  is 
extrapolated  with  zero  gradient  for  the  inlet.  At  the  exit,  constant  pressure  is  set  while  the  velocities 
are  extrapolated  with  zero  gradient.  It  should  be  noted  that  the  inflow  boundary  condition  only 
holds  for  the  initial  processor  used  in  the  calculation.  Similarly,  the  exit  boundary  condition  applies 
only  to  the  last  processor  assigned  in  the  parallel  processing.  Additional  boundary  conditions  on 
both  sides  of  each  subdomains  are  a  function  of  time  and  neighboring  domains. 


Numerical  Results 

A  total  of  five  3-D  unsteady  calculations  are  presented.  The  Space  Shuttle  Main  Engine’s  (SSME) 
coaxial  injector  has  been  chosen  for  baseline  modeling  case.  First,  the  injector  of  the  SSME  main 
combustion  chamber  (MCC)  has  been  investigated.  Further  simulations  have  been  performed  for 
the  MCC  injector  with  a  longer  recess  length  and  lower  velocity  ratio  in  order  to  assess  the  impact 
of  these  parameters.  In  addition,  a  simulation  of  SSME  prebumer  injector  has  been  conducted. 
The  both  injectors  use  liquid  oxygen  (LOX)  and  gaseous  hydrogen  ( GH2 )  as  an  oxidizer  and  fuel 
respectively,  but  have  different  geometries  and  physical  properties  consistent  with  the  fuel-rich 
prebumer  gases  entering  the  main  combustion  chamber.  Figure  4  illustrates  the  schematic  of  three 
dimensional  coaxial  injector  with  geometries  of  main  burner  and  prebumer  injectors.  The  physical 
properties  of  each  injector  flow  in  operating  condition  are  provided  in  Table  1. 

SSME  main  burner  coaxial  injector 

The  nondimensionalization  was  based  on  the  geometries  in  Figure  4  and  Table  1  by  using  fuel 
annulus  diameter  Df,  liquid-phase  velocity  Ui  and  liquid  density  p;  as  nondimensional  parameters. 
Under  this  nondimensionalization,  the  LOX  post  inner  diameter  has  a  value  of  0.55,  the  post  tip 
thickness  (t^ox)  is  0.05  and  the  recess  length  ( L )  is  0.75.  The  present  result  assumes  a  gas/liquid 
density  ratio  of  0.02  consistent  with  high  pressure  combustion  conditions.  As  discussed  in  the  grid 
refinement  study,  a  mesh  size  of  72  x  71  x  71  is  used  for  this  case.  With  parallel  processing, 
12  processors  are  used  and  the  computational  mesh  is  split  into  12  sets,  resulting  8  x  71  x  71 
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of  subdomain  mesh  assigned  to  each  processor.  Typically,  a  150,000  time  step  calculation  with  a 
timestep  At  =  0.0001  (0.3msec)  takes  about  three  weeks. 

The  results  exhibit  characteristics  termed  ”self-pulsation”  and  ”self-oscillation”  as  theorized  by 
Bazarov  16”17.  Pulsations  are  evident  in  changes  in  the  flowrate  with  time,  while  oscillations  are 
evidenced  by  azimuthal  motion  of  the  central  liquid  core  about  the  annulus.  Prior  2-D  simulations 
19)20  have  shown  similar  behavior,  although  the  azimuthal  motion  cannot  be  resolved  with  this 
decreased  degree  of  freedom. 

The  overall  density  field  behavior  of  the  liquid  jet  for  conditions  roughly  equivalent  to  the 
SSME  MCC  injector  is  illustrated  in  Figure  5.  Unfortunately,  the  model  could  not  be  run  at  the 
gas/liquid  density  ratio  of  this  injector,  but  converged  calculations  were  demonstrated  at  a  density 
ratio  approximately  double  that  in  the  actual  engine.  Other  simulations  at  the  lower  density  ratio 
were  obtained  and  will  be  discussed  subsequent  to  this  case.  The  left  column  in  Fig.  6  depicts 
density  contours  in  a  cutaway  view  interpretation  of  the  motion  of  the  jet.  The  right  column  shows 
the  density  contours  at  the  exit  plane. 

The  resultant  highly  nonlinear,  quasi-periodic  oscillation  occurs  naturally  as  a  result  of  the 
Kelvin-Helmholtz  instability  mechanism.  There  is  no  oscillation  in  the  inflow,  yet  strong  oscil¬ 
lations/pulsations  develop  as  a  result  of  the  large  relative  dynamic  pressure  between  the  gas  and 
liquid  streams.  The  development  of  self-pulsation  mode  of  the  liquid  jet  is  apparent  in  those  fig¬ 
ures. 

Pressure  contour  plots  at  times  corresponding  to  those  in  the  density  plots  are  shown  in  Figure 
6.  Initially,  we  start  the  calculations  assuming  uniform  pressure  in  the  recessed  region  assessed  in 
the  computations  as  shown  in  the  t=0  contour.  Nonlinear  wave  growth  occurs  quite  rapidly  and  a 
highly  3-D  field  develops  after  initial  trasient  period.  High  pressure  regions  develope  in  accordance 
with  azimuthal  and  axial  motion  of  the  structures  on  the  jet.  Longitudinal  wave  motion  is  apparent 
in  the  low  pressures  at  times  t=8.15,  8.75  vs.  the  high  pressure  at  t=8.45.  Multiple  high  pressure 
regions  can  be  formed  at  different  azimuthal  locations  within  the  gas  annulus  as  evidenced  in  the 
t=8.90  and  t=14.60  images. 

Velocity  streamlines  of  the  jet  in  x  —  y,  x  —  z  and  y  —  z  planes  are  shown  in  Figure  7  along 
with  velocity  contours  at  exit  plane.  Only  three  time  frames  corresponding  to  the  last  three  in  each 
density  plot  are  shown  here.  The  asymmetric  flow  recirculations  in  the  region  just  downstream  of 
the  LOX  post  are  shown  in  left  and  middle  columns.  The  far  right  column  illustrates  axial  velocity 
contours  at  the  exit  plane.  While  vortices  shed  from  the  separated  flow  region  at  the  base  of  the 
LOX  post  are  evident,  the  velocity  contours  show  that  the  jet  comes  out  of  the  orfice  exit  at  an 
average  velocity  across  the  entire  exit  plane  except  near  the  wall  and  flow  recirculation  area.  This 
agrees  with  the  experiment  of  Mayer,  Schik  and  Schaffler 6  which  verified  the  rapid  vanishing  of 
relative  velocity  between  the  central  liquid  jet  and  annular  gaseous  jet  at  the  downstream  of  LOX 
post  tip  area.  In  the  right  column  of  the  figure,  streamlines  at  the  orifice  exit  plane  are  depicted  to 
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indicate  the  complex  cross-currents  which  develop  in  the  channel. 

One  interesting  bulk  measure  of  orifice  performance  is  the  massflow  delivered  at  the  orifice 
exit  plane  as  a  function  of  time.  This  quantity  can  be  computed  by  quadrature  of  the  density- 
axial  velocity  product  over  the  exit  plane  area.  Figure  8  shows  the  time  history  of  mass  flow 
for  the  density  field  computed  in  Fig.  5  and  the  attendant  velocity  field.  The  unsteadiness  (self¬ 
pulsations)  of  the  jet  are  quite  evident  in  this  figure;  massflow  variations  of  39%  about  the  mean 
flow  are  apparent. 

Spectral  analysis  of  the  data  is  performed  in  a  manner  common  for  experimentaldata.  Because 
the  massflow  rate  was  recorded  only  for  every  hundredth  iteration,  aliasing  is  expected  to  be  a 
problem.  That  is,  just  as  in  experiments,  there  is  reason  to  believe  that  physical  processes  are 
happening  at  frequencies  in  excess  of  the  Nyquist  frequency  for  the  data  stream.  Thus,  an  anti¬ 
aliasing  filter  was  convolved  with  the  signal,  the  the  signal  was  windowed,  and  the  mean  subtracted 
out.  The  power  spectrum  was  computed  through  use  of  an  FFT.  Spectra  computed  without  this 
method  showed  significantly  worse  signal  to  noise  ratio  when  compared  to  those  presented  here. 
The  spectral  content  of  the  mass  fluctuation  signal  provides  information  regarding  the  frequency 
spectrum  present.  Using  a  Fast  Fourier  Transform,  the  spectral  content  is  shown  in  Figure  9. 
There  is  a  significant  energy  content  over  a  range  of  dimensionless  frequencies  with  a  maximum 
at  /  =  2.2074,  which  represent  around  73,000  Hz  in  dimensional  units.  This  frequency  is  about 
triple  the  frequency  at  which  liquid  is  replaced  in  the  submergence  region. 

Along  with  the  primary  frequency  at  /  =  2.2074,  the  first  harmonic  frequency  at  /  =  4.4816 
is  also  shown  in  this  figure.  This  higher  frequency  corresponds  to  the  small  peak  variation  at  the 
top  of  primary  mass  flow  fluctuation  signal  in  Figure  8,  which  is  caused  by  small  surface  vortices 
coming  out  of  the  nozzle  exit  in  an  irregular  manner. 

In  a  recent  study  of  Smith  and  Mayer  6,  the  liquid  core  jet  oscillation  in  a  coaxial  injector  was 
observed,  and  the  frequency  of  flame  flickering  in  a  combustion  chamber  was  reported  as  near  7000 
Hz,  which  falls  within  the  same  frequency  band  of  the  SSME  main  burner  coaxial  injector’s  case. 
The  same  velocity  and  density  ratios  were  used  with  similar  geometry.  Even  though  their  coaxial 
injector  was  not  exactly  the  same  one  used  in  this  study,  their  experimental  study  provides  some 
evidence  that  the  liquid  jet  oscillation/pulsation  plays  a  significant  role  in  combustion  response 
within  the  chamber.  This  is  an  important  finding  requiring  further  experimental  verification. 

SSME  MMC  injector:  lower  velocity  ratio 

A  case  of  lower  velocity  ratio,  Ug/Ui  =  6,  has  been  investigated  in  order  to  assess  the  influence  of 
this  parameter.  As  the  gas-phase  velocity  is  reduced  to  180m/s,  the  corresponding  Reynolds  num¬ 
ber  of  the  gas-phase  is  4.5  x  105,  but  the  density  ratio  of  default  value  (0.01)  for  SSME’s  MCC  was 
used  for  this  simulation.  The  time  history  of  mass  flow  is  plotted  in  Figure  10.  A  small  magnitude 
self-pulsation  appears  during  the  initial  time  period  but  it  decays  gradually  without  showing  any 
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surface  vortices  or  wave-type  instability.  The  spectral  content  of  the  damping  fluctuation  signal  is 
depicted  in  Figure  1 1.  A  primary  frequency  at  /  =  2.2704  with  very  small  amplitude  is  observed 
in  this  case;  a  value  similar  to  that  computed  for  the  baseline  case. 

In  prior  2-D  simulations  19,20 ,  reductions  in  the  relative  gas/liquid  velocity  also  showed  a  sta¬ 
bilizing  influence  since  the  dynamic  pressure  forcing  the  instability  scales  with  the  square  of  this 
velocity.  However,  the  conclusion  that  the  Ug/Ui  =  6  case  is  stable  for  this  injector  is  somewhat 
tenuous  as  our  present  simulation  assumes  a  completely  steady  flow  entering  both  streams.  A 
more  comprehensive  analysis  or  confirmation  via  experimental  data  would  be  required  to  confirm 
the  relative  stability  of  this  (or  any)  case. 

SSME  MCC  injector:  longer  recess  length 

The  recess  length  (depth)  is  an  important  paramater  governing  this  flow  process,  so  a  calculation 
was  performed  to  assess  the  impact  of  lengthening  the  recess  by  50%  ( L  =  1.25)  keeping  all  other 
parameters  fixed  as  in  the  baseline  MCC  calculation  described  previously.  In  order  to  keep  the 
same  grid  aspect  ratio  and  to  stay  within  a  reasonable  run-time,  a  mesh  of  132  x  75  x  75  is  used 
with  22  processors.  With  the  increased  number  of  processors,  each  processor  handles  the  same 
mesh  size  of  subdomain  (8  x  75  x  75)  as  set  in  the  baseline  case.  This  methodology  provided  an 
acceptable  run-time  of  about  three  weeks. 

Density  contours  in  3-D  cylinder  and  at  exit  plane  are  depicted  in  Fig.  12.  The  development 
of  self-pulsations  is  very  clear  from  the  beginning  and  it  continues  to  be  strong  through  the  entire 
modeling.  Azimuthal  oscillations  become  apparent  in  latter  part  of  simulation  as  shown  in  exit 
plane  density  contours  in  the  left  column.  The  liquid  core  moves  more  actively  in  an  irregular 
pattern  than  the  baseline  case.  In  this  simulation,  the  liquid  jet  shows  a  more  structured  wave-type 
motion  in  the  axial  direction.  The  additional  flowpath  length  permits  the  longer  wave  instability 
to  become  more  pronounced;  the  internal  liquid  jet  is  even  seen  to  break  up  within  the  submerged 
region  as  shown  in  image  (b)  in  the  left  hand  column  of  the  figure. 

The  last  three  time  frames  in  density  contour  plots  are  used  to  show  velocity  streamlines  and 
axial  velocity  contours  at  the  orifice  exit  plane  in  Figure  13.  Flow  recirculations  are  shown  in  x  —  y 
and  x  -  z  plane  in  the  first  and  second  columns  while  velocity  contours  at  exit  plane  are  depicted  in 
the  third  column  with  velocity  streamlines  in  cross-stream  directions.  Vortex  structures  shed  from 
the  base  of  the  LOX  post  are  also  evident  in  this  case. 

Figures  14  and  15  show  the  time  history  of  mass  flow  and  the  spectral  content  of  this  signal, 
respectively.  The  magnitude  of  pulsation  is  larger,  showing  oscillating  components  of  the  order 
of  48%  of  the  mean  flow.  The  frequency  of  oscillation  is  lower  than  that  observed  in  the  baseline 
geometry.  The  larger  magnitude  pulsation  can  be  explained  by  the  development  of  wave-type 
oscillation  due  to  the  longer  recess  length.  The  primary  frequency  is  decreased  to  /  =  1.1352,  but 
this  value  is  still  consistent  with  the  prior  observation  in  that  it  is  still  roughly  triple  the  frequency 
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at  which  the  liquid  traverses  the  submerged  region.  Similar  trends  have  been  observed  in  prior 
2-D  simulations  19>2°;  increasing  submergence  length  increases  the  time  available  for  instabilities 
to  grow. 

SSME  pre-burner  coaxial  injector 

The  SSME’s  prebumer  geometry  is  shown  in  Figure  4  and  flow  properties  are  tabulated  in  Table 
1.  The  prebumer  injector  has  a  shorter  recess  length  (L  =  0.5)  and  three  times  thicker  LOX  post 
thickness  ( tLox  =0.15)  than  main  burner  injector.  The  LOX  post  diameter  (Dt)  is  0.45,  reducing 
the  mixture  ratio  of  liquid  to  gas-phase  significantly.  The  density  ratio  is  much  higher  (e  =  0.04) 
due  to  the  higher  pressures  in  this  combustor,  but  the  velocity  ratio  is  identical  to  the  SSME  MCC 
baseline  case.  The  same  mesh  size  of  72  x  71  x  71  is  used  for  the  prebumer  case,  so  the  axial 
resolution  is  somewhat  better  due  to  the  reduced  submergence  length. 

The  simulation  showed  signs  of  weak  pulsation/oscillation  activity  during  the  initial  transient, 
but  these  oscillations  were  damped  and  an  essentially  undisturbed  jet  issued  from  the  orifice  at  long 
times.  In  prior  2-D  simulations  19-21 ,  increasing  gas  density  magnified  the  aerodynamic  interac¬ 
tion,  resulting  in  a  faster  and  finer  atomization.  However,  thickening  the  LOX  post  and  reducing 
the  submergence  length  both  had  stabilizing  influences  on  the  instability.  In  this  case,  those  sta¬ 
bilizing  effects  overwhelmed  the  instability  caused  by  the  aerodynamic  interaction  at  higher  gas 
density,  therefore  the  overall  jet  appears  as  stable. 

As  noted  previously,  the  present  simulations  assume  idealized  conditions  at  the  entry  to  the 
submerged  region  and  that  the  real  device  may  indeed  exhibit  unsteadiness.  More  comprehensive 
simulations  would  be  necessary  to  fully  illuminate  this  issue. 

Main  burner  injector  geometry  with  preburner  flow  properties 

Having  the  previous  result  with  prebumer  injector  case  raises  questions  about  relation  between 
flow  properties  and  injector  geometry.  In  order  to  investigate  the  reciprocal  effects  between  those 
parameters,  a  simulation  has  been  conducted  using  the  MCC  injector  geometry  with  the  prebumer 
flow  properties.  In  summary,  the  recess  length  L  is  0.75  with  LOX  post  diameter,  D{  =  0.55  and 
thickness,  tLOx  =  0.05.  The  density  ratio  e  is  0.04  and  velocity  ratio  stays  at  12.  The  numbers  in 
Table  1  are  used  for  the  Reynolds  numbers  for  gas  and  liquid-phase.  The  Reynolds  numbers  for  this 
case  should  have  been  changed  based  on  the  main  burner  injector  geometry  with  the  prebumer  flow 
properties,  but  the  same  values  of  the  prebumer  case  has  been  used  in  error.  However,  the  Reynolds 
numbers  are  in  the  same  order  of  magnitude  even  with  the  change  and  it  was  decided  that  the  time 
invested  in  the  expensive  calculation  was  such  that  the  values  used  would  be  suitable.  In  prior  2-D 
parametric  studies  19>2°,  Reynolds  number  had  only  a  minor  influence  on  results.  Therefore,  the 
result  of  this  case  is  still  valid  to  predict  qualitative  behavior  of  the  liquid  jet. 


145 


Figure  16  illustrates  the  limit-cycle  behavior  of  the  density  contours  in  the  passage.  Pulsations 
and  oscillations  similar  to  those  in  the  baseline  MCC  simulation  are  observed.  One  interesting 
observation  from  this  case  is  that  the  pulsations/oscillations  become  weaker  in  magnitude  after  a 
certain  transient  period.  This  behavior  is  evidenced  in  the  orifice  massflow  history  shown  in  Fig. 
17.  The  main  difference  between  this  case  and  the  MCC  baseline  case  discussed  in  Figs.  5-9  is 
that  the  density  ratio  is  twice  as  large.  In  prior  2-D  simulations,  increased  gas  density  served  to 
increase  the  magnatude  and  frequency  of  the  instability  presumably  due  to  the  fact  that  the  gas 
can  impart  more  inertia  to  the  liquid.  Here,  in  the  3-D  case,  we  see  a  reduction  in  pulsations  with 
increased  gas  density. 

The  spectral  content  of  the  waveform  shown  in  Fig.  18  are  shown  in  Figure  19  provides  a 
partial  explanation  for  the  reduction  in  magnatude  of  the  oscillation  for  this  higher  gas  density 
case.  Significant  energy  content  is  stored  in  a  number  of  frequencies;  i.e.  the  oscillation  is  not 
organized  at  a  single  frequency  so  its  overall  amplitude  is  lower  than  the  baseline  MCC  case.  The 
same  primary  frequency  that  was  observed  in  the  baseline  main  burner  case  is  present  (/  =  2.0), 
but  in  general  the  fluctuation  is  very  irregular  and  unpredictable.  Further  parametric  studies  are 
required  to  illuminate  this  complex  issue. 


Summary  and  Conclusions 

A  3-D  incompressible,  unsteady,  viscous  Navier-Stokes  solver  for  two-phase  flow  has  been  devel¬ 
oped  for  full  three  dimensional  simulations  of  coaxial  injectors  in  liquid  rocket  engines.  The  Space 
Shuttle  Main  Engine’s  (SSME)  main  burner  and  prebumer  injectors  served  as  baseline  modeling 
cases. 

High  amplitude  hydrodynamic  instabilities  are  observed  in  the  recessed  region  where  com¬ 
putations  are  made.  The  resulting  frequencies  and  massflow  variations  for  the  5  simulations  are 
summarized  in  Table  2.  The  parameter  whose  effect  was  investigated  is  underlined.  The  massflow 
of  the  device  oscillates  as  much  as  48%  about  the  mean  flow  for  the  cases  studied  indicating  a 
strong  presence  of  ’’self  pulsations”  as  theorized  by  other  researchers.  Oscillations  of  the  liquid 
core  jet  about  the  gas  annulus  are  also  observed  in  some  of  the  cases.  Reduction  in  the  length  of 
the  recess  region  or  the  velocity  difference  between  the  streams  serves  as  a  stabilizing  influence. 
Increases  of  gas  density  introduced  complex  effects  in  terms  of  spectral  content  and  amplitude  of 
the  massflow  oscillations. 

For  the  SSME  main  combustion  chamber  injector,  the  primary  harmonic  frequency  was  about 
7300  Hz.  Recent  experiments  at  DLR  at  comparable  conditions  indicated  flickering  of  the  flame  at 
a  similar  frequency  of  7000  Hz.  For  this  reason,  there  may  be  strong  connections  between  hydro- 
dynamic  instabilities  and  the  resultant  combustion  processes.  Additional  studies  are  recommended 
for  complete  verification  of  this  important  conclusion. 
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Nomenclature 

D  -  Channel  width 
L  -  Channel  length 

h  -  Liquid  sheet  thickness  from  center  line 
twx  -  Liquid  post  thickness 

A U  -  Velocity  difference  between  liquid  and  gas  phase 

e  -  Density  ratio  ( pg/pi ) 

m  -  Nozzle  mass  flow  rate 

P  -  Pressure 

Re  -  Reynolds  number 

a  -  Void  fraction 

fi  -  Viscosity 

p  -  Fluid  pseudo-density 

CPi  -  Computational  processor 


Subscripts 

i  -  Inlet 
o  -  Outlet 
l  -  Liquid-phase 
g  -  Gas-phase 
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Figure  Captions 


1.  Schematic  of  domain  splitting  used  in  parallel  processing. 

2.  Cross-sectional  view  of  3-D  computational  mesh:  Axial  direction. 

3.  Grid  refinement  study  on  three  different  mesh  sizes  (SSME  main  burner  condition:  e  =  0.02, 
Ug/Ui  =  12,  L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105) 

4.  Schematic  of  SSME  coaxial  injector  and  summary  of  critical  dimensions. 

5.  Cutaway  view  (left  column)  and  orifice  exit  plane  (right  column)  of  density  contours  showing 
typical  Self-Pulsation:  SSME  Mainbumer  injector,  e  =  0.02,  Ug/Ut  =  12,  L/D  =  0.75, 
Rei  =  1.1  x  106,  Reg  =  9.0  x  105 

6.  Cutaway  view  (left  column)  and  orifice  exit  plane  (right  column)  of  pressure  contours  in  3-D 
cylinder  at  various  times:  SSME  Mainbumer  injector,  e  =  0.02,  Ug/Ut  =  12,  L/D  =  0.75, 
Ret  =  1.1  x  106,  Reg  =  9.0  x  105 

7.  Velocity  streamlines  in  cross-sectional  view  (left  and  middle  columns)  and  orifice  exit  plane 
(right  column):  SSME  Mainbumer  injector,  e  =  0.02,  Ug/Ut  =  12,  L/D  =  0.75,  Ret  = 
1.1  x  106,  Reg  =  9.0  x  105 

8.  Time  history  of  mass  flow  at  exit  plane:  SSME  Mainbumer  injector,  e  =  0.02,  UgjUx  =  12, 
L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105 

9.  Spectral  analysis  of  mass  flow  fluctuation:  SSME  Mainbumer  injector,  e  =  0.02,  Ug/Ut  = 
12,  L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  9.0  x  10s 

10.  Time  history  of  mass  flow  at  exit  plane:  lower  velocity  ratio,  e  =  0.01,  Ug/Ui  =  Q,L/D  = 
0.75,  Ret  =  1.1  x  106,  Reg  =  4.5  x  105 

11.  Spectral  analysis  of  mass  flow  fluctuation:  lower  velocity  ratio,  e  =  0.01,  Ug/Ut  =  6, 
L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  4.5  x  105 

12.  Density  field  evolution  for  injector  with  a  longer  recess  length:  e  =  0.01,  Ug/Ut  =  12, 
L/D  =  1.25,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105 

13.  Velocity  streamlines  and  contours  with  a  longer  recess  length:  e  =  0.01,  Ug/Ut  =  12, 
L/D  =  1.25,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105 

14.  Time  history  of  mass  flow  at  exit  plane:  longer  recess  length,  e  =  0.01,  Ug/Ut  =  12, 
L/D  =  1.25,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105 
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15.  Spectral  analysis  of  mass  flow  fluctuation:  longer  recess  length,  e  =  0.01,  Ug/Ui  =  12, 
L/D  =  1.25,  Rei  =  1.1  x  106,  Reg  =  9.0  x  105 

16.  Density  contours  in  initial  development  of  Self-pulsation:  e  =  0.04,  Ug/Ui  =  12,  L/D  = 
0.75,  Ret  =  6.6  x  105,  Reg  =  3.37  x  106  post  thickness,  tLOx  =  0.1 

17.  Density  contours  in  small  magnitude  Self-pulsation  mode:  e  =  0.04,  Ug/Ui  =  12,  L/D  = 
0.75,  Rei  =  6.6  x  105,  Reg  =  3.37  x  106 

18.  Time  history  of  mass  flow  at  exit  plane:  main  burner  geometry  with  prebumer  flow  proper¬ 
ties,  e  =  0.04,  Ug/Ui  =  12,  L/D  =  0.75,  Ret  =  6.6  x  105,  Reg  =  3.37  x  106 

19.  Spectral  analysis  of  mass  flow  fluctuation:  e  =  0-04,  Ug/Ui  =  12,  L/D  =  0.75,  Rei  = 
6.6  x  105,  Reg  =  3.37  x  106 


r 


152 


Table  1:  Physical  properties  of  SSME  coaxial  injectors  in  operating  conditions 


Physical  properties 

SSME  mainbumer 

SSME  prebumer 

Chamber  pressure  (MPa) 

19.3 

34.1 

Liquid  injection  velocity  (m/s) 

31.3 

30 

Liquid  density  ( kg/m 3) 

1117 

1125 

Liquid  Reynolds  No. 

1.1  x  106 

6.60  x  105 

Gas  injection  velocity  (m/s) 

360.6 

360 

Gas  density  (kg/m3) 

9.47 

42.5 

Gas  Reynolds  No. 

9.0  x  105 

3.37  x  106 

Density  ratio  (pg/pi) 

0.0085 

Velocity  ratio  ( Ug/Ui ) 

12 

12 

Mixture  ratio  (liq./gas) 

0.89 

0.98 
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Table  2:  Summary  of  3-D  simulations  results 


CASE 

e 

UJU, 

L/D 

tLOX 

A  m/m 

frequency 

MCC 

0.02 

12 

0.75 

0.05 

39 

2.2  (7300Hz) 

MCC-U 

0.01 

6 

0.75 

0.05 

Near  0 

2.27  (7500Hz) 

MCC-L 

0.01 

12 

1.25 

0.05 

48 

1.14  (3800Hz) 

PB 

0.04 

12 

05 

0.15 

Near  0 

N/A 

12 

0.75 

0.05 

24 

0.5/1. 4/2.0/4.7/5.1 
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Grid  Refinemnet  Study 


Figure  3:  Grid  refinement  study  on  three  different  mesh  sizes  (SSME  main  burner  condition:  e 
0.02,  Ug/Ui  =  12,  L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105 
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Injector  Geometric  Parameters 

SSME  Prebumer 

SSME  Mainbumer 

Fuel  Annulus  Diameter,  Df 

5.03mm 

8.84mm 

LOX  Post  Outer  Diameter, D0 

3.76mm 

6.60mm 

LOX  Post  Inner  Diameter,  Dj 

2.26mm 

4.77mm 

Post  Tip  Thickness,  tLOX 

0.76mm 

0.92mm 

Recess  Length,  L 

2.54mm 

6.48mm 

Figure  4:  Schematic  of  SSME  coaxial  injector  and  summary  of  critical  dimensions 


Figure  6:  Cutaway  view  (left  column)  and  orifice  exit  plane  (right  column)  of  pressure  contours  in 
3-D  cylinder  at  various  times:  SSME  Mainbumer  injector,  e  =  0.02,  Ug/Ui  =  12,  L/D  =  0.75, 
Rei  =  1.1  x  106,  Re„  —  9.0  x  10s 


Figure  10:  Time  history  of  mass  flow  at  exit  plane:  lower  velocity  ratio,  e  =  0.01,  Ug/Ui  =  6, 
L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  4.5  x  105 


Figure  11:  Spectral  analysis  of  mass  flow  fluctuation:  lower  velocity  ratio,  e  =  0.01,  Ug/Ui  —  6, 
L/D  =  0.75,  Ret  =  1.1  x  106,  Reg  =  4.5  x  105 


Density  Contours  in  3-D  cylinder  &  Exit  plane 


Figure  12:  Density  field  evolution  for  injector  with  a  longer  recess  length:  e  =  0.01,  u,/u, 
L/D  =  1.25,  Rei  =  1.1  x  106,  Reg  =  9.0  x  105 


0.8 


0.7 


time 


Figure  14:  Time  history  of  mass  flow  at  exit  plane:  longer  recess  length,  e  =  0.01,  Ug/Ui  =  12, 
L/D  =  1.25,  Ret  =  1.1  x  106,  Reg  =  9.0  x  105 


Figure  15:  Spectral  analysis  of  mass  flow  fluctuation:  longer  recess  length,  e  =  0.01,  Ug/Ui  =  12, 
L/D  =  1.25,  Rei  =  1.1  x  106,  Reg  =  9.0  x  105 
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ity  contours  in  initial 
=  6.6  x  105,  Reg  = 


Figure  17:  Density  contours  in  small  magnitude  Self-pulsation  mode:  e  =  0.04,  Ug/Ui 
L/D  =  0.75,  Ret  =  6.6  x  105,  Reg  =  3.37  x  106 


Figure  18:  Time  history  of  mass  flow  at  exit  plane:  main  burner  geometry  with  prebumer  flow 
properties,  e  =  0.04,  Ug/Ut  =  12,  L/D  =  0.75,  Ret  .=  6.6  x  10s,  Reg  =  3.37  x  106 


Figure  19:  Spectral  analysis  of  mass  flow  fluctuation:  e  =  0.04,  Ug/Ui  =  12,  L/D  =  0.75, 
Rei  =  6.6  x  105,  Reg  =  3.37  x  106 


168 


